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1.0 INTRODUCTION 


Global warming and population growth have resulted in an increase in 
the intensity of natural (e.g., climate change) and anthropogenic stressors. 
Investigating the responses of environmental processes to the cumulative 
effects of stressors — now known as cumulative effects assessment (CEA) 
- was already being practiced (e.g., groundwater resources described in 
Schoff & Sayan, 1969) when the concept of environmental impact assess- 
ment (EIA) was introduced in 1969. However, the theory and concept of 
CEA was defined by the US Council on Environmental Quality (CEQ) 
in 1978. The concept of CEA has then been gradually described in more 
detail by other scholars (e.g., Canter, 1999; Ross, 1998; Cooper, 2004), and 
many well-designed approaches have been proposed in the literature (e.g., 
Dubé & Munkittrick, 2001; Dubé et al., 2013; Lokke et al., 2010). 

There are different definitions for CEA in the literature (Cooper & 
Sheate, 2002; Bérubé, 2007; Noble et al., 2014). However, it is defined here 
as an assessment of cumulative environmental changes due to human and 
natural stressors over a period of time for the past, present, and future, rela- 
tive to a baseline or standard. CEA not only includes analyzing and mod- 
elling environmental changes, but it also supports planning alternatives 
that promote environmental monitoring and management. A variety of 
methods are used in literature (Table 1) for cumulative effects assessment, 
such as spatial analysis, network analysis, interactive matrices, and eco- 
logical modelling (Smit & Spaling, 1995). In addition to demanding the 
development of models to help answer policy questions (Castronova et al., 
2013), CEA also requires observations (i.e., monitoring) of changes in nat- 
ural phenomena. 

Investigating the complex nature of environmental problems re- 
quires the integration of different environmental processes across major 


components of the environment, such as water, climate, ecology, air, and 
social aspects. The increasing dissatisfaction resulting from disjointed and 
narrowly focused environmental management approaches has recently 
encouraged the use of integrated environmental modelling approaches 
(Jakeman & Letcher, 2003). Integrated environmental modelling refers to 
coupling of thematic based numerical or conceptual models to solve com- 
plex real-world problems involving the environment and its relationship 
to human systems and activities. 

The concept of developing integrated models first appeared decades 
ago (Mackay, 1991). However, there has recently been an increased inter- 
est in further developing integrated modelling frameworks in response 
to the emergence of problems related to regional scale land-use manage- 
ment, impacts of global climate change, evaluation of ecosystem services, 
fate and transport of nanomaterials, and life-cycle analysis. A variety of 
models have been developed to investigate the processes of each individ- 
ual environmental component and the way they interact with each other. 
However, they have failed to consider environmental processes of other 
components of the environment and their complex interplay within the 
environment as a whole. Integrated modelling frameworks are often the 
only way to take into account the important environmental processes, 
interactions, relevant spatial and temporal scales, and feedback mechan- 
isms of complex systems for CEA. The other obstacle is the uncertainty as 
to whether an applied modelling system meets its intended purposes and 
sufficiently represents reality, an issue which is reinforced by a lack of clear 
understanding of models’ mechanisms. In this regard, this book looks at 
(i) understanding interactions and relationships between environmental 
components, such as climate, land, hydro, ecology, and resulting responses 
due to anthropogenic/natural stressors in CEA, (ii) reviewing modelling 
approaches for each component, (iii) reviewing existing integrated model- 
ling systems for CEA, and (iv) proposing an integrated modelling frame- 
work and perspectives for future research avenues for CEA (Figure. 1). 
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Figure 1. Overview of the Components of this Report 
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1.0 Introduction 


Table 1. Examples of Cumulative Effects Studies 


Authors Objective Study area Method Indicators Response 
Shrestha et Quantify Athabasca Soil and Blue Results pro- 
al., 2017 the impacts River Basin Water (surface jected the 
of climate (ARB), Assessment freshwater climate of 
change on Canada Tool & ground- the ARB to 
monthly, (SWAT), water) and be wetter by 
seasonal, calibrated green water 21-34% and 
and annual over 1990- (soil mois- warmer by 
water 2005 period ture) 2.0-5.4 °C 
balances on an annual 
of blue and time scale. 
green water The annual 
resources at average blue 
sub-basin, and green 
regional, water flow 
and basin- (streamflow 
wide spatial and evapo- 
scales. transpira- 
tion) was 
projected 
to increase 
by 16-54% 


and 11-34%, 
respectively. 


Cho et al., Simulate Athabasca Community Annual Average 

2017 dry, wet River Basin Multi-Scale dry, wet, nitrogen 
and total (ARB), Air Quality and total deposition 
deposition Canada (CMAQ) depositions increases 
of acidifying model of acidifying rom the 
compounds sulphur and historical to 
using four nitrogen existing and 
different compounds uture cases. 
emission Sulphur 
scenarios deposition 
and provide decreases 
guidance rom the 
on possible historical 
priorities o existing 
for oil sands cases but 
emissions increases to 
manage- uture cases 
ment for even though 
addressing regional SO, 
acid depo- emission 
sition in the continuously 
Alberta’s oil decreases. 
sand region 
(AOSR). 
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Table 1. (continued) 


Authors Objective Study area Method Indicators Response 
Zhang etal., Provide a Baker Time series Flow regime The 
2016 quantitative Creek and cross- magnitude, 
assess- Willow River correlation variability, 
ment on watersheds, analysis and and return 
how forest British paired-year period of 
disturbances Columbia, approach high flows 
affect the Canada in the Baker 
components Creek 
of flow watershed 
regimes at a were 
large water- increased 
shed scale. on average 
by 154.3%, 
324.2%, and 
11 years, 
respectively, 
and the 
timing 
of high 
flows was 
advanced by 
about 9 days 
during the 
disturbed 
periods. 
Ahmed, Discern and Rideau River Numerical Low and It was 
2013 quantify the watershed, modelling peak demon- 
hydrologic Ontario, techniques discharge strated that 
functions Canada (Mike 11) the flood 
of wetlands risk would 
within the increase if 
context of wetlands are 
the Rideau removed. 
River The low flow 
watershed. will likely 
increase if 
wetlands are 
removed. 
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Table 1. (continued) 


Authors Objective Study area Method Indicators Response 


Deitch etal., Predict Sonoma GlIS-based Streamflow Results 
2013 streamflow County, hydrologic illustrate 
impairment California, model that 
caused by USA impairment 
438 small caused by 
reservoirs reservoirs 
in the study varies 
area appreciably 
over space, 
but as 
reservoirs fill 
over time, 
impairment 
is lower 
through 
most of the 
drainage 
network. 


Cormier et Assessing Big Darby Conceptual Fish and A decrease 
al., 2000 ecological Creek models for molluscs, in the 
riskina watershed, exposure water quali- percentage 
watershed central Ohio, of fish and ty, sedimen- of Tanytar- 
USA benthic tation, and sini midges 
macroin- hydrologic along with 
vertebrates, regimes an increase 
surveys of in the per- 
fish and cent of tox- 
molluscs, ic-tolerant 
biological invertebrate 
indices species 
found in fish 
is expect- 
ed to be 
consistently 
associated 
with an in- 
crease inthe 
concentra- 
tions of toxic 
chemicals in 
the water. 
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2.0 INTERACTIONS BETWEEN 
ENVIRONMENTAL COMPONENTS FOR 
CUMULATIVE EFFECTS ASSESSMENT 


2.1 Climate and hydrology 


The hydrological cycle is linked with the Earth’s radiation balance (Figure 
2). Changes in radiation balance due to increasing amounts of greenhouse 
gases tend to increase atmospheric and oceanic temperatures and alter 
the seasonal frequency, intensity, and location of precipitation. Changes 
in precipitation patterns can cause a change in runoff and consequently 
influence soil moisture and infiltration. 

In addition, an increase in temperature can result in a change in air 
humidity and alter the evaporation and evapotranspiration processes. 
The melting of snow and glacier ice that results in rising global sea levels 
is also caused by the increased temperature that would lead to increased 
winter runoff and decreased late spring and summer runoff. An increase 
in the evapotranspiration rate will decrease soil water storage and will 
change the amount of water that infiltrates and percolates into the soil. 
Climate change can also alter the type of native vegetation cover, the leaf 
area index (LAI), and the soil moisture content at local and regional scales 
(Arnell, 1994). 

Temperature, precipitation, and evaporation are the main climate 
variables that impact water resources at local to global scales and are 
common climate inputs for hydrological models. Temperature is the most 
important variable compared to other meteorological components, since 
along with its effects on hydrological processes, it has a significant impact 


on precipitation and evaporation. An increase in temperature increases 
rain in proportion to snow in areas where precipitation is a combination 
of snow and rain in the winter season. 

An increase in precipitation generally increases the amount of water 
stored in the soil and thus has an influence on recharge to groundwater, 
actual evapotranspiration, and surface runoff. Surface runoff supplies riv- 
ers and lakes; this water, along with rainfall, will then recharge ground- 
water. Spatial and temporal changes in precipitation will likely alter tim- 
ing, duration, and magnitude of aquifer recharge, thus affecting ground- 
water in return. Changes in flood frequency are influenced by the timing 
and duration of rainfall (and/or snowmelt) events, and changes in drought 
frequency may occur as a result of a lack of rainfall. 

Evaporation is directly influenced by temperature and precipitation. 
It increases with an increase in temperature, since evaporation occurs in 
response to the availability of thermal energy. A decrease in precipitation 
results in increased solar radiation (leading to more frequent occurrences 
of clear skies). Consequently, there is a decrease in humidity and therefore 
increased evaporation. Evaporation is thus considered a key link between 
the atmosphere and the soil water within the hydrological cycle. 

Studies that examine the relationship between climate and hydrology 
began in the early 1940s (Veijalainen, 2012) but were limited due to avail- 
ability of computing technology. In the 1980s, with drastic advances in the 
domain of climate science and hydrological modelling, a major step was 
taken in the study of climate systems and their interaction with hydrol- 
ogy. However, there was still considerable uncertainty in the calculations 
due to unreliable scenarios and low model resolution (Veijalainen, 2012). 
In the 1990s, many comprehensive and complex climate and hydrological 
models were developed, and as computers became more powerful, studies 
began to probe into different aspects of climate and hydrology. 
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Figure 2. Impacts of Climate Change on Hydrological Processes 
(modified from Arnell, 1994) 


frre sana na sn nn nn nn nn nnn nnn nnn n nanny 


Hydrological response to climate 


Climate Change 
[Greenhouse Gas (GHG)] 


Increase in absorption 
and emitting of radiation 


Groundwater storage -- Surface water storage 


2.0 Interactions Between Environmental Components 


Recently, several studies have been conducted to evaluate climate change 
impacts on hydrology (Boyer et al., 2010; Cherkauer & Sinha, 2010; 
Grillakis et al., 2011; Forbes et al., 2011; Stefan et al., 2012; Barron et al., 
2012; Farjad et al., 2016; Bajracharya et al., 2018). In general, most of these 
studies have considered the following approach (Di Baldassarre et al., 
2011): 


1. use of a scenario or set of global climate model (GCM) 
scenarios; 


2. transfer the scale of projected GCMs variables into local 
scale (downscaling); 


3. apply the downscaled GCMs data into hydrological 
models; and 


4. simulate responses of hydrological processes to climate 
change. 


Different studies have been carried out in each of the above steps of the 
framework using different climate and hydrological models (Matondo et 
al., 2004; Steele-Dunn et al., 2008; Kuhn et al., 2011; Tshimanga & Hughes, 
2012; Dobler et al., 2012). Each study focused on different aspects such as 
hydrological processes of interest, required climate variables, geography 
of the study area, scale, and availability of models and data. These studies 
have illustrated changes in hydrological processes in response to climate 
and could be grouped in three broad categories: 


1. Changes in runoff/streamflow: A number of studies have 
explored potential changes in mean annual runoff in 
response to climate change (Jones et al., 2006; Gardner, 
2009). Chang and Jung (2010) assessed the effects of 
climate change on annual, seasonal, and high and low 
runoff in the 218 sub-basins of the Willamette River 
basin in Oregon. The study showed that snowmelt- 
dominated basins exhibit large reductions in summer 
flow in response to increased temperature, while rainfall- 
dominated basins show large increases in winter flow 
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in response to precipitation change. Thus, an increase 
in temperature and precipitation can affect the quantity 
and seasonal distribution of streamflow (Campbell et 
al., 2011). Jha et al. (2004) found that a unit increase in 
precipitation due to climate change will cause a larger 
increase in streamflow. 


Changes in subsurface water (including soil water, 
deeper-vadose zone water, and confined and unconfined 
aquifer waters): Any impact of climate change on surface 
hydrology has a potential influence on the subsurface 
hydro-geology as both systems are coupled together and 
inseparable. However, feedback between soil water and 
unsaturated zones with climate change occurs over short 
periods of time compared to saturated zones, which 
might take decades (Green et al., 2011). Goderniaux et 

al. (2009) estimated the potential climate impacts on 
groundwater in the Geer Basin in eastern Belgium. They 
applied a physically based surface-subsurface flow model 
combined with six climate change scenarios. They found 
that the groundwater level is expected to decrease up to 
8 m by 2080. 


Change in extreme events: A change in climate will 
result in a change in the magnitude of events and their 
frequency of occurrence. The extreme hydrology- 
related events can be analyzed in two categories. The 
first is peak flow (flood) frequency, which refers to the 
peak discharges for recurrence intervals. Mareuil et 

al. (2007) assessed flood frequency and severity in the 
Chateauguay River Basin, in the province of Quebec, 
Canada, and they recognized that both traits are affected 
by spring snowmelt and summer/fall storms, since an 
increase in temperature will have considerable effect on 
earlier spring runoff in areas where snowmelt composes 
the major proportion of the streamflow. The second 
category is low flow (drought) frequency, which refers 
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to the lowest streamflow during dry periods for a given 
recurrence interval. Yulianti and Burn (1998) investigated 
the influence of air temperature on low flow frequency 
due to climate change in the Prairies Region of Canada. 
They found that the frequency of low flow events has 
increased due to an increasing temperature trend over the 
past century, and that there is a significant relationship 
between regional air temperature and low flows in the 
Prairies. 


2.2 Land, ecology, and climate 


Changes in land surface/ecology could have effects on regional and global 
climate variables (Li & Molders, 2008; Steyaert & Knox, 2008; Findell 
et al., 2009). Many studies have shown that land practices such as de- 
forestation (Meher-Homji, 1991; Werth & Avissar, 2002) and agriculture 
(Feddema et al., 2005) may affect regional climate to a similar or greater 
extent than would climate change driven by global changes in atmospher- 
ic chemistry alone. However, land-related variations in surface exchanges 
can have a more pronounced impact at a local scale than globally. In their 
study, Findell et al. (2007) concluded that observed changes in land sur- 
face properties have slight impact on globally averaged climatic variables, 
while such changes are highly significant locally in the annual mean and 
in most months of the year. 

Changes in land surface cause changes in the surface energy and 
moisture balances (Findell et al., 2007) and consequently can impact cli- 
mate variables such as air temperature, precipitation, and evapotranspir- 
ation. Also, Werth and Avissar (2002) found that deforestation can result 
in a significant reduction of summertime precipitation in the Amazon 
region as well as the surrounding area. Changes in the physical character 
of the land surface include surface roughness, soil texture, and vegetation. 
Changes in surface properties will result in changes in albedo (the propor- 
tion of the incident light or radiation that is reflected by a surface), leafarea 
index (LAI, leaf area per unit ground surface area), and root depth (Figure 
3). Different surface properties have different influences on climate. For 
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instance, forests have more net radiation (the balance between incoming 
and outgoing energy that can be exchanged between land surface and cli- 
mate) than pastures. This is because pastures have a higher albedo with 
more longwave radiation but a lower aerodynamic roughness, as well as 
less evapotranspiration (shallow-rooted vegetation cannot extract water 
in deep soil) compared to forests not only during the dry season, but also 
in the wet season due to the reduced roughness and available energy. 
Therefore, forests can make the climate cooler and pastures can make it 
warmer. 


Figure 3. Changes in the Physical Character of the Land Surface 
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Even though evapotranspiration in forests generally exceeds that of grass- 
land based on an annual time scale, this might not be true in all climate 
conditions and temporal scales due to the complex behavior of guard cells 
around stomata. It is likely that the strong regulation of stomatal open- 
ing in response to radiation, temperature, vapour pressure deficit, and 
the larger rooting depth contributes to the conservative character and 
persistence of evapotranspiration in forests. Therefore, it is possible that 
evaporative cooling over grasslands exceeds that over forests when there is 
plenty of soil moisture; however, the opposite could happen when the soil 
moisture condition is low (Teuling et al., 2010). 

The relative division of available energy between sensible and latent 
heat fluxes (types of energy released or absorbed in the atmosphere) is 
sensitive to land cover conditions, as is the division of precipitation be- 
tween runoff and evaporation (Pitman, 2003; Findell et al., 2007). This 
is because less latent heat results in less water vapour to the atmosphere 
and tends towards decreasing cloudiness and precipitation. Decreases in 
sensible heat flux tend to cool the planetary boundary layer (Figure 4) and 
lessen convection (Betts et al., 1996; Pitman, 2003). 

Variations in the vegetation cover change the surface area of vegeta- 
tion that is in contact with the atmosphere and the balance between fluxes 
from the soil and vegetation (Pitman, 2003). Changes in the leaf area index 
(LAI) can result in variations in the exchange of both sensible heat flux 
and latent heat flux (Figure 5). LAI can also affect the interactions between 
the biosphere and the atmosphere through various processes such as res- 
piration, transpiration, photosynthesis, and rain interception. 

The root zone has the ability to return soil water to the atmosphere 
through evapotranspiration. Root depth (RD) can affect the extent of the 
exchange between soil and the atmosphere. A change in climate can re- 
sult in changes in vegetation and consequently the root depth, which can 
ultimately affect the global water cycle. RD is adversely affected by the 
removal of natural vegetation, which reduces moisture availability leading 
to increased sensible heat flux and warmer and deeper boundary layers. 
Conversely, reforestation usually increases the latent heat flux via tran- 
spiration and interception loss (Betts et al., 1996; Pitman, 2003; Findell et 
al., 2008). The amount of soil moisture available to plants for transpira- 
tion can undergo dramatic variations if the distribution of roots changes. 
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Furthermore, one may also find positive feedback between reduced root 
water uptake, rainfall, and further reductions in root depth (Figure 6). 

Changes in climate variables can influence the vegetation species 
and vegetation phenology. As mentioned, extremely cold or hot soil tem- 
peratures can damage plant growth. Warmer climates can cause earlier 
flowering and accelerate the growth of plants. Extremely cold climates can 
damage cell walls of growing plants in spring or fall. On the other hand, 
drier climates and lack of precipitation can damage plant cells and roots, 
which could lead to leaf scorch, wilting, and eventually leaf drop. Wetter 
climates with too much water can reduce the amount of oxygen in the soil, 
which can result in root loss as well as vulnerability to fungal diseases. 
Climate change and variability can also influence animal species; for ex- 
ample, an increase in sea temperature can result in a loss of oxygen, which 
can negatively affect certain fish species. 


Figure 4. Responses of Vegetation and Climate Components to a 
Change in Albedo (modified from Pitman, 2003) 
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Figure 5. Responses of the Hydrological Cycle to a Change in LAI 
(modified from Pitman, 2003) 


G 
Change in 
LAI 
Change in canopy Change in stomatal Change in 
shading conductance interception loss 


[ 


V ° 

Change in surface Change in Change in 3 
energy absorption transpiration throughfall x 
4 vo 

an 

4 ; 

6 


----------> 


Change in Change in soil > Change in 
infiltration moisture evaporation 


16 


2.0 Interactions Between Environmental Components 


Figure 6. Responses of the Hydrological Cycle to a Change in RD 
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2.3 Land use/land cover (LULC) and hydrology 


In addition to climate variables, changes in land use/land cover (LULC) 
can also impact hydrological processes. LULC changes alter water bal- 
ance, causing changes in infiltration, evapotranspiration, and intercep- 
tion. LULC changes result from anthropogenic activities such as urbaniz- 
ation, agriculture, deforestation, and afforestation. These activities modify 
the physical properties of the surface such as vegetation and soil, which in 
turn alter impervious surfaces, soil moisture, and runoff. This ultimately 
leads to changes in surface and groundwater balance and may cause floods 
and/or a change in drought frequency (Arnell, 1994). For example, the 
conversion of natural vegetation into agricultural land decreases evapo- 
transpiration and hence increases fresh water availability (Gordon et al., 
2003), while the conversion of land and natural vegetation into impervi- 
ous surfaces decreases infiltration, evapotranspiration, and groundwater 
recharge and increases the rate of surface runoff (Arnold & Gibbons, 
1996). Arnold and Gibbons (1996) reported that the conversion of forest 
into 10-20%, 35-50%, and 75-100% impervious surfaces increases runoff 
two times, three times, and five times, respectively. Consequently, urban- 
ization changes many of the processes that affect streamflow by replacing 
vegetation and soil with impervious surfaces (McGrane, 2016). Rose and 
Peters (2001) identified five major impacts of urbanization on hydrology: 
(i) a higher proportion of precipitation converts to surface runoff while 
impervious surfaces increase; (ii) the catchment response to precipita- 
tion is accelerated and the lag time between precipitation and runoff is 
decreased; (iii) peak flow magnitudes are increased for all but the largest 
storm events; (iv) low flow is decreased due to reduced contributions 
from groundwater storage; and (v) water quality is degraded by effluent 
discharges. 

The impacts of LULC changes on hydrological processes can be sum- 
marized as follows: 


Runoff and infiltration: Many studies have demonstrated 
the relationships between LULC properties and runoff and 
infiltration. Different results have been obtained in different 
studies, but most of them demonstrated an inverse relation 


18 2.0 Interactions Between Environmental Components 


between forest cover and runoff for many types of forested 
landscapes (S. Wang et al., 2008). Bosch and Hewlett (1982) 
concluded that changing annual and perennial herbaceous 
vegetated land by establishing forest cover decreases water 
yield. In this order, coniferous hardwood, brush, and grass 
cover have the highest consumptive use of water with 

the implication that conversions from these respective 
cover types to annual crops will increase water yield 
accordingly. Bruijnzeel (1990) concluded that vegetation 
type, plant maturation stage, and management practices 
affect the change in water yield. Calder (1998) reported that 
vegetation after clear-cutting, depending on the species, 
may consume more water than very old forests. Other 
studies indicate deforestation increases runoff (Coe et al., 
2011), water yield (Bosch & Hewlett, 1982), and flood peaks 
(Hornbeck et al., 1970; Harr et al., 1975; Harr, 1981, 1986). 


Baseflow: LULC changes are one of the main human- 
induced activities affecting the groundwater system. 

LULC changes that lead to increased infiltration, such as 
terracing, contour farming, and organic matter build-up, 
may offset the adverse effects caused by deforestation. 

For example, agricultural practices can offset the adverse 
effects of deforestation by increasing infiltration resulting 
in increased baseflow (Bruijnzeel, 1990). Conversion of land 
for agricultural purposes alters key vegetation parameters 
that affect recharge, including fractional vegetation 
coverage, wilting point, and rooting depth. Reducing 
fractional vegetation coverage to zero during fallow periods 
between crop rotations can increase recharge as shown in 
the Northern Great Plains of the United States (Scanlon et 
al., 2005). 


Streamflow/Peak flows: Peak flow is significantly inversely 
related to the infiltration and water storage capacity of soil. 
Verry et al. (1983) and Bruijnzeel (1990) found increase 
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in peak flows after cutting trees in a watershed. Also, 
Bruijnzeel (1990) reported the diminishing role of soil 

and plants in modulating peak flows resulting from larger 
magnitude of storm precipitation. In large watersheds, the 
effects of land-use changes on peak flows may not be that 
pronounced at the watershed outlet mainly due to the time 
lag between different tributaries and spatial variation in 
rainfall (Bruijnzeel, 1990); however, it is more significant 
in smaller sub-watersheds (Brooks et al., 2003). According 
to Calder (1998), increasing transpiration in dry periods 
will increase soil moisture deficits and reduce dry season 
flows. It is widely reported that an increase of urban lands 
is usually associated with an increase in high streamflow, 
decrease in low streamflow, and an increased variability in 
streamflow. This is a result of the increase in impervious 
surface caused by urbanization, which decreases infiltration 
of precipitation and increases runoff (White & Greer, 2006; 
Tu, 2009). 


2.4 Air quality, climate, water, and ecology 


Air pollutants are emitted from a range of both natural (e.g., wind-blown 
dust, wildfires, volcanos) and manmade sources including mobile sources 
(e.g., transportation vehicles), stationary sources (e.g., power plants, in- 
dustrial facilities, and factories), and aerial sources (e.g., cities, agricul- 
tural area, and wood burning fireplaces). The six common air pollutants, 
known as the criteria pollutants, are particulate matter (PM), nitrogen 
dioxide (NO,), sulfur dioxide (SO,), carbon monoxide (CO), ozone (O,), 
and lead (Pb). Air pollutants can either be released directly into the atmos- 
phere as primary emissions or formed as a result of chemical interactions 
between precursor substances. CO, Pb, NO,, and SO, are directly emitted 
from different air pollution sources. PM can be formed when emissions 
of SO,, NO, ammonia, organic compounds, and other gases react in the 
atmosphere, or it can be emitted. Ozone is not directly emitted but is 
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formed when oxides of nitrogen (NO,) and volatile organic compounds 
(VOCs) react in the presence of sunlight. 

Gases and particulates from the atmosphere can fall into water, soil, 
and on vegetation in dry or wet forms by a process known as air depos- 
ition or atmospheric deposition. Dry deposition is a slow and continuous 
process, whereas wet deposition can occur at a faster rate when precipita- 
tion intercepts gases and particulates and falls on the landscape. However, 
both processes depend on several other environmental factors and are 
influenced by the properties of the gases or particles (Table 2). 

Air pollution, especially that caused by emissions of sulphur and 
nitrogen and ground-level ozone, can influence ecosystems. Sulphur diox- 
ide and nitrogen oxides can deposit as acid rain on the landscape (e.g., on 
water bodies, soil, and vegetation). This can result in increasing water and 
soil acidity with adverse effects on the ecosystem. For example, soil acid- 
ity can alter the soil chemistry, which can consequently influence water 
quality and plant growth. Atmospheric deposition of nitrogen can be leth- 
al for aquatic organisms such as fish. It is a risk for the eutrophication (the 
process of accumulation of nutrients) of water systems. Nutrient overload 
in aquatic ecosystems can cause algae blooms (a rapid accumulation in 
the population of algae in marine water bodies or freshwater) and con- 
sequently a loss of oxygen. An increase in ground-level ozone can result in 
the loss of plant cover by damaging plant cell membranes, inhibiting key 
processes. 


Table 2. Factors Affecting Dry and Wet Deposition of Gases and 
Particles 


Components Wet deposition Dry deposition 

Aerosol particles Cloud parameters Particle properties 
Precipitation Near-surface concentration 
Particle properties Weather condition 

Gases Gas-specific parameters Soil, surface, vegetation properties 
Cloud parameters Physical and chemical properties 
Precipitation of tracer 


Near-surface concentration 
Chemical reactions 
Weather condition 
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3.0 MODELLING APPROACHES FOR 
EACH COMPONENT 


3.1 Hydrological models 


Hydrological models are used to quantitatively analyze hydrological pro- 
cesses. They are applied to understand the dynamic interactions between 
the climate and land surfaces by simulating complex interaction processes 
in the hydrological cycle, which are subject to both normal and extreme 
climate conditions (Zhang, 2007) or to changes in the physical character- 
istics of a land surface. 

The study of hydrological modelling dates back to the 1850s when 
Mulvany (1850) used a rational method to calculate the volume of runoff 
based on the percentage of rainfall and watershed area (Zhang, 2007). 
Sherman (1932) developed a unit hydrograph concept in which the runoff 
process was assumed to be linear and time invariant (Dawdy, 1983). The 
first runoff model based on physical processes was developed by Horton 
(1933) using the infiltration-excess runoff theory (Zhang, 2007). In the 
1950s, Kalinin and Milyukov (1958) developed a channel routing method 
using linear analysis (Dawdy, 1983). 

The classification of models helps to understand their capability and 
structure, while each class is an individual representation of the hydro- 
logical cycle. There is no universal model to simulate and characterize 
watershed hydrology; models are classified in different ways depending 
on the criteria of interest. Selection of these models is based on the spatial 
and temporal scale of the studies, the type of watershed, the modelling 
objective (Zhang, 2007), data availability, and economic constraints. 
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This current study distinguishes models that treat a watershed as a 
spatially variable system in three classes: lumped, semi-distributed, and 
distributed models. 


3.1.1 Lumped models 

Lumped models consider the watershed as one computational unit, with 
state variables representing average values of watershed characteristics, 
such as the total rainfall, soil moisture, or overland flow. The derivation 
of such variables depends on empirical relationships derived by various 
techniques, including curve fitting, using available monitoring data. 

Lumped models rely on the techniques of systems analysis in relating 
inputs to outputs without reference to the internal and physical mechan- 
isms of the watershed. Calibration of the model is based on the comparison 
between observed and simulated watershed outflows. There are different 
lumped models with different functional forms defined intuitively. These 
models are commonly called conceptual models in hydrology. 

Generally, flow-routing mechanisms over the watershed area are ig- 
nored in lumped models (Beven, 2001). The discharge of lumped models 
is based on the global dynamics of the system. These models ignore in- 
filtration of surface runoff and its connection with river flow since they 
are not physically based. They require many assumptions that increase 
the uncertainty of the models. For instance, precipitation is considered 
uniformly distributed over the watershed spatially and temporally. LULC, 
soil, and geology are also assumed uniform across a watershed (Reed, 
2004). Although they have some advantages such as having a simple struc- 
ture and easy setup, calibration, and use, they require long-term historical 
data for calibration and the parameter values may be potentially difficult 
to physically interpret. 

Numerous studies have used different lumped models for hydrologic- 
al studies. Four examples of lumped models are listed in Table 3 and de- 
scribed as follows: 
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IHACRES: Croke and Jakeman (2007) applied the 
IHACRES model that was originally designed for temperate 
climates to assess streamflow in three different climatic 
regions in Australia. They concluded that the model is 
suitable for arid and semi-arid catchments. However, they 
also suggested that the length of the calibration period 
should be increased to accommodate the lower frequency of 
streamflow events. 


SRM: Kustas et al. (1994) used three different approaches 
for modelling snowmelt: 1) a degree-day model called the 
snowmelt runoff model (SRM); 2) a restricted degree-day 
model, characterized by a simple radiation component 
combined with the degree-day approach to improve 
estimates of snowmelt and reduce the need to adjust the 
melt factor over the ablation season; and 3) a daily energy 
balance model. They tested the three approaches using melt 
rates with lysimeter outflow measurements. The restricted 
degree-day and energy balance models produced better 
results than the snowmelt runoff model. 


WATBAL: Yates (1996) investigated the impact of climate 
change on discharge in Arkansas. He selected two different 
basins to evaluate the range of the model’s applicability; one 
(Mulberry River) in a humid climate that was dominated 
by winter rainfall and warm summers and a second one 
(East River) in a semi-arid region that was dominated by 
snowfall and colder temperatures. WATBAL as a water 
balance model combined with the Priestley-Taylor method 
were used to estimate potential evapotranspiration. The 
WATBAL model was originally designed as a simple model 
to assess the impact of climate change on a watershed. The 
parameters of the model are direct runoff, surface runoff, 
subsurface runoff, maximum catchment water-holding 
capacity, and baseflow. The results revealed that the model 
behaves fairly well, given its simplicity. The model showed 
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the sensitivity to precipitation change in the Mulberry 
Basin. They suggested that WATBAL lacks seasonal 
parameters, as there was a strong seasonal variation in 
runoff in the Mulberry Basin. 


USDAHL: England (1975) investigated soil moisture in 
two layers of soil in an Oklahoma basin. He simulated soil 
moisture for two arbitrary layers, 0 to 9 inches and 10 to 
33 inches, and then compared the modelled results with 
observed data during a 15-month period. The results of 
the model simulation of soil moisture were very close to 
the observed data in layer 1, but there was a large deviation 
between simulated results and observed data in winter and 
spring for layer 2. He thereby concluded that the model is 
capable of simulating soil moisture continuously at a site. 
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Table 3. Lumped Models 


Model Provider Reference Description Input data Output 
IHACREC Center for Jakeman The model Rainfall, Streamflow, 
Ecology and et al., employs unit temperature, wetness 
Hydrology 1990 hydrograph evapotrans- index 
(CEH) (UH) and piraton 
Wallingford, simulates 
UK, and the steamflow 
Australian either 
National continuously 
University or individually 
(ANU) 
SRM Swiss Martinec, The model Precipitation, Streamflow 
Snow and 1975 (simple temperature, 
Avalanche degree-day)is and elevation 
Research designed for 
Institute basins where 
(SSARI) snowmelt is a 
major runoff 
component 
WATBAL University Yates, The model is Precipitation, PET, 
Helsinki, 1994 a soil water temperature, discharge 
Finland balance humidity, 
model. sunshine 
duration 
USDAHL United States Holtan et The model Precipitation, streamflow, 
Department al., 1975 is designed temperature, AET, soil 
of Agriculture to simulate evapotranspi- moisture, 
Hydrograph continuous ration groundwater 
Laboratory streamflow recharge 


predictions. 
It was useful 
to evaluate 
interactions 
between 
agricultural 
activities and 
hydrology of 
small rural 
watersheds. 
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3.1.2 Semi-distributed models 

Semi-distributed models consider conceptual functional relationships for 
homogeneous sub-catchments as lumped units. These models discretize 
landscape based on common land use, soil, and slope characteristics of 
a watershed, known as hydrologic response units (HRUs). They include 
some of the important features of a watershed compared to the lumped 
models and require less data, and they have lower computational costs 
compared to distributed models (Orellana et al., 2008). However, HRUs 
are often spatially disconnected and routed directly to sub-basin outlets. 
Table 4 lists examples of semi-distributed models. 


3.1.3 Distributed models 
Distributed models consider catchments as finite geo-referenced com- 
putational units with different responses to forcing inputs. A grid-based 
hydrological model may not necessarily be a distributed model, unless 
grid cells can interact both vertically and horizontally with adjacent cells 
(within surface, unsaturated, or saturated zones) to simulate surface 
and/or subsurface processes. The main category of distributed models is 
physically based models, which are defined in terms of theoretically con- 
tinuum equations based on physics. Although distributed models require 
large amounts of data for parameterization, they provide more detail of 
hydrological processes (Refsgaard, 1996) and imply a discrete grid system 
in which the spatial variations are aggregated over each grid. These mod- 
els are used for: 1) flood studies - evaluation of volume and timing of peak 
flows; 2) yield studies - evaluation of the total flow obtained from rainfall 
within a watershed; 3) low flow studies — assessment of low flow in the 
watershed; and 4) water management studies - assessment of the impacts 
of the manmade water management infrastructure. 

Ina different classification, Refsgaard (1996) recognized four practical 
applications for using distributed physically based models: 


1. Watershed changes: Distributed physically-based models 
are suited to estimating the influence of both natural 
and manmade changes on the hydrological cycle, since 
the parameters of the models have direct physical 
interpretation. Predictions by these models are based 
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Table 4. Semi-Distributed Models 


Model Provider Reference Description Input data Output 
HBV-96 Swedish Lindström et The model Precipita- Daily 
Meteoro- al., 1997 is designed tion, tem- flow, soil 
logical and for runoff perature, moisture 
Hydrologic simulation and monthly ET 
Institute forecasting. 
(SHMI) 
SLURP National Kite, 1997 The model is a Precipita- Infiltration, 
Hydrological daily timestep tion, tem- overland 
Research hydrological perature, flow, soil 
Institute, model dividing albedo, moisture 
Saskatoon, a watershed snowmelt 
Canada into a number rate 
of units known 
as Aggregated 
Simulation 
Areas (ASAs). 
HEC-HMS US Army Kumar & The model is Precipita- Infiltration, 
Corps of Bhattacharjya, designed to tion, tem- soil 
Engineers 2011 simulate both perature, moisture, 
(US-ACE) event and evapotrans- runoff, flow 
Hydrologic continuous piration in open 
Engineering simulation over channels 
Center long periods of 
time. 
SWAT USDA Preksedis et The model Precipita- Surface 
Agricultural al., 2008 predicts the tion, tem- and 
Research effect of perature, subsurface 
Service management soil types runoff, flow 
at the decisions routing 
Grassland, on water, through 
Soil and sediment, drainage 
Water nutrient and network, 
Research pesticide yields ET, soil 
Laboratory with reasonable moisture, 
in Temple, accuracy on 
TX, USA large, ungauged 
river basins. 
VIC University of Liang et al., Variable Meteo- Runoff, soil 
Washington 1994; Wood Infiltration rological moisture, 
et al., 2004; Capacity (VIC) parameters ET, and 
Zhao etal., is a semi- (e.g., pre- infiltration 
2013; Lievens distributed, cipitation, 
et al., 2016 grid-based air tem- 
macroscale perature, 
model designed and wind 
for large-scale speed) and 
hydrological soil 
modelling. 
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on parameter values and physical characteristics of a 
watershed, such as land use and soil. 


Simulations with intensive and short-term records: In 
contrast to lumped models, which require long historical 
data for the assessment of the parameters, distributed 
models simulate hydrological processes in the watershed 
using short records. 


Prediction of the ungauged watershed: In a well-gauged 
watershed, the model can be calibrated against observed 
discharge data with less uncertainty, but a prediction of 
flows in an ungauged watershed is complex, with a higher 
level of uncertainty regarding model output. The physical 
significance of the model parameters allows distributed 
physically-based models to predict the responses of an 
ungauged watershed. 


Spatial assessment: Distributed models can use spatially 
variable inputs and predict outputs spatially. They provide 
necessary information such as movements of rainstorms, 
groundwater abstractions, and recharge, while lumped 
models can only consider average values in a watershed. 


Table 5 lists four examples of distributed models identified in the liter- 
ature. Some of the hydrological models have been used in combination 
with water quality modelling or have been integrated with groundwater 
models. A more detailed description of these hydrological model applica- 
tions appears below: 
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WATFLOOD: Toth et al. (2006) applied WATFLOOD model 
to investigate the hydrological regimes of the Peace and 
Athabasca watersheds under five different climate change 
scenarios (GCMs). Results revealed a significant shift 
towards an earlier melt season, a shift in the timing of peak 
flows, and small changes in the annual flow volumes. 
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CASC2D: Marsik and Waylen (2006) applied the two- 
dimensional, physically-based hydrologic model CASC2D 
to evaluate the influence of land-use/land cover change on 
hydrology from 1979 to 1999 in the Quebrada Estero in 
Costa Rica. The results showed increased peak discharges 


and above-threshold flood durations with changing LU/LC. 


This model was found to be well suited for operational use 
in tropical watersheds like the Quebrada Estero. 


MIKE SHE/MIKE 11: Farjad et al. (2016) investigated 
seasonal and annual responses of hydrological processes 
to climate change in the 2020s and 2050s in the Elbow 
River watershed, southern Alberta, Canada. The MIKE 
SHE/MIKE 11 model was applied to simulate different 
hydrological processes under the GCM-scenarios of 
NCARPCM-A1B, CGCM2-B2(3), HadCM3-A2(a), 
CCSRNIES-AIFI, HadCM3-B2(b)). The model was set up 
based on a rigorous sensitivity analysis along with three 
different methods of calibration and validation to capture 
the complex watershed hydrology. Results indicated that 
future climate change is expected to progressively modify 
hydrological processes over the next 60 years. 


HydroGeoSphere: Davison et al. (2018) assessed streamflow 
characteristics in the downstream and upstream of 

the Athabasca River Basin in Alberta. They found that 
forestlands and peatlands have a strong influence on the 
hydrology of the watershed. 
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Table 5. Distributed Models 


Model Provider Reference Description Input data Output 
CASC2D US Army Re- Julien et Fully Precip- Soil moisture, 
search Office al., 1995 unsteady, itation, infiltration, 
(ARO) funded two-dimen- tempera- surface and 
Center for sional, infiltra- ture, ET channel 
Excellence in tion-excess runoff 
Geoscience at (Hortonian) 
Colorado hydrologic 
State Univer- model 
sity 
WATFLOOD University Kouwen, The model Precip- Surface 
of Waterloo, 2001 is designed itation, low, soil 
Canada for real-time tempera- moisture, ET, 
flood fore- ture, ET nfiltration 
casting 
MIKE SHE/ DHI Wijesekara The model is Precip- Overland 
MIKE 11 et.al.,2012 anintegrated _ itation, lows, base 
hydrological tempera- lows, AET, 
modelling ture, ET soil moisture, 
system for nfiltration, 
building and groundwater 
simulating able 
surface water 
flow and 
groundwater 
flow 
HydroGeo- Aquanty Hwang et It is an Precip- Streamflow, 
Sphere al., 2018 integrated itation, AET, 
surface water tempera- Infiltration, 
and ground- ture, ET groundwater 
water model table 
to simu- 
late major 
hydrological 
processes 


3.2 Water quality models 


Water quality is a complex subject and may involve interactions between 
surface water, groundwater, and coastal water systems. It is controlled by 
characteristics of watersheds and aquifers, including climate, land cov- 
er and land uses, geology, lithology, chemical reactions, and anthropo- 
genic activities (Praskievicz & Chang, 2009; Tsakiris & Alexakis, 2012; 
Whitehead et al., 2009). Ever since Streeter and Phelps (1925) built the well- 
known Streeter-Phelps Oxygen Sag Formula (SP model) to describe the 
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oxygen balance in the Ohio River, the development of water quality mod- 
els has undergone tremendous improvements. Many water quality models 
have been developed to tackle various water quality problems or threats 
associated with growing populations, urbanization, and industrialization. 
Models have evolved from a primary stage (prior to 1965), characterized 
by simple BOD-DO bilinear systems considering point source pollution, 
to an improving stage (1965-1995), when two- and three-dimensional 
nonlinear system models were built to include non-point source pollution 
inputs, capacities of hydrodynamic, sediment, eutrophication simulation, 
and linkages to watershed models (Q. Wang et al., 2013). Thereafter, in- 
tegrated modelling systems at various levels of sophistication have been 
developed where atmospheric deposition, climate and land use changes, 
and interactions between surface water, groundwater, and water resource 
management are considered in modelling to evaluate their impacts on wa- 
ter quality (Burian et al., 2002; Hesse & Krysanova, 2016; Hien et al., 2015; 
Panagopoulos et al., 2015; Wellen, Kamran-Disfani, & Arhonditsis, 2015). 

Water quality models are effective tools to simulate and predict the 
transport and fate of pollutants in aquatic environments and support en- 
vironmental impact assessment and planning. They can be classified ac- 
cording to their characteristics and intended purposes (Sharma & Kansal, 
2013), for example: 


e modelling purpose (simulation, optimization), 
e development (generic, site specific), 
e model type (physical, mathematical), 


e application area (rivers, lakes, reservoirs, watershed, 
groundwater, estuaries, integrated), 


e constituents of concern (sediments, salts, nutrients, 
metals, PAHs, etc.), 


e nature (deterministic, stochastic), 
e spatial variation (1-, 2-, 3-dimensional), 


e spatial resolution (lumped, semi-distributed, 
distributed), 
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e temporal variation (steady state, quasi-dynamic, 
dynamic simulation), or 


e solution method (analytical, finite difference, finite 
element, linear, nonlinear and dynamic programming, 
etc.). 


There are already many attempts to review the development and applica- 
tions of receiving water quality models (Bahadur, Amstutz, & Samuels, 
2013; Tsakiris & Alexakis, 2012) and watershed water quality models 
(Booty & Benoy, 2009; K. H. Cho et al., 2016; Wellen et al., 2015). While 
most reviews are presented in brief articles for selected models, ASCE 
(2017) published a book offering a comprehensive review of 13 watershed 
models and 13 receiving water quality models, for total maximum daily 
load (TMDL) development. Building on these previous efforts, this book 
summarizes the advances in water quality modelling techniques with a 
focus on integrated modelling. It has been recognized that peer-reviewed 
journal literature no longer provides a representative picture of the subject 
of regional integrated environmental modelling, as modelling systems are 
becoming “too big to be published” or too pragmatic (Barthel & Banzhaf, 
2016; Wood, 2012). Therefore, the existing reviews and water quality mod- 
elling studies reported in journal publications, conference proceedings, 
technical reports, and software manuals since 2000 have been searched 
and screened according to relevance to watershed and integrated water 
quality modelling. In the following subsections, water quality modelling 
approaches, along with the characteristics and applications of the models, 
are described. 


3.2.1 Receiving Water Quality Models 
The receiving water quality models commonly also include hydrodynamic 
models. Models vary widely as to their water quality modelling capacities, 
as they are developed with various temporal and spatial variables, types 
of receiving waters, contaminants of interest, and representations of pro- 
cesses (ASCE 2017). 

Review of selected receiving water quality models are summarized in 
Table 6 according to key model characteristics or features, followed by a 
brief description about the source, capabilities, and applicability for each. 
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CE-QUAL-W2 

CE-QUAL-W2 is a two-dimensional, longitudinal/vertical, hydrodynam- 
ic, and water quality model for stratified and non-stratified rivers, estuaries, 
lakes, reservoirs, and river basin systems. The model has been under con- 
tinuous development since 1975. The original model was known as LARM 
(Laterally Averaged Reservoir Model) developed by Edinger and Buchak 
(1975). The first LARM application was on a reservoir with no branches. 
Subsequent modifications that allowed for multiple branches and estuarine 
boundary conditions resulted in the code known as GLVHT (Generalized 
Longitudinal-Vertical Hydrodynamics and Transport Model). Addition 
of the water quality algorithms by the Water Quality Modelling Group at 
the United States Army Engineer Waterways Experiment Station (WES) 
resulted in CE-QUAL-W2 version 1.0 (Environmental and Hydraulics 
Laboratory, 1986). The latest release is version 4.2.2, released in August 
2020 and distributed by Portland State University. The model and source 
code are publicly available at http://www.ce.pdx.edu/w2/. 

The CE-QUAL-W2 software directly links a hydrodynamic module 
and a water quality module using a dynamic coupling approach. In the 
model, the geometry of a waterbody is represented by a finite difference 
computation grid defined using layers of segments and cells. The module 
predicts water surface elevations, velocities (longitudinal and vertical), 
and temperatures. Temperature is included in the hydrodynamic calcula- 
tions because of its effect on water density. The model can calculate onset, 
growth, and breakup of ice cover. Water quality computations are done 
after a hydrodynamic computation, allowing for feedback between water 
quality and hydrodynamic variables. The effects of salinity or total dis- 
solved solids on density and thus on hydrodynamics is simulated only if it 
is set up as one of the state variables of the water quality module. The water 
quality algorithm is modular, allowing constituents to be easily added as 
additional subroutines. The model simulates eutrophication, alkalinity, 
and generic water quality processes in water column and sediment diagen- 
esis processes in the sediment bed. Water quality variables include generic 
constituents, sediments, nutrients, multiple algal groups, epiphyton, peri- 
phyton, zooplankton, macrophytes, carbonaceous biochemical oxygen 
demand, dissolved oxygen, and dissolved and particulate labile/refractory 
organic matters. The generic water quality groups can be used to define 
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any number of conservative tracers, water age or hydraulic residence time, 
coliform bacteria, and contaminants. Additionally, more than 60 derived 
variables, such as pH, TOC, DOC, TON, TOP, DOP, TP, TN, TKN, and 
turbidity, can be computed internally from the state variables and the out- 
put can be compared to measured data. 

As the water surface elevation is solved implicitly, eliminating the 
surface gravity wave restriction on timestep, CE-QUAL-W2 permits lar- 
ger timesteps during a simulation, resulting in decreased computational 
time. As a result, the model can easily simulate long-term water quality 
responses (Cole & Wells, 2017). Note that water quality can be updated less 
frequently than hydrodynamics, thus reducing computational require- 
ments. However, water quality is not decoupled from the hydrodynam- 
ics (i.e., separate, standalone code for hydrodynamics and water quality), 
as output from the hydrodynamic model is stored on disc and then used 
to specify advective fluxes for the water quality computations. The CE- 
QUAL-W2 model is a powerful and widely used laterally averaged longi- 
tudinal/vertical two-dimensional model for simulating hydrodynamics 
and water quality in rivers, lakes, reservoirs, and estuaries (ASCE, 2017; 
Shabani, Zhang, & Ell, 2017). The model has been further enhanced to de- 
velop the Cumulative Environmental Management Association (CEMA) 
Oil Sands Pit Lake Model. This development incorporated a sediment 
diagenesis module, tailings consolidation, pore water release, biogenic gas 
production, bubble release, and salt rejection during ice formation with- 
in the CE-QUAL-W2 model Version 3.6 (Berger & Wells, 2014; Prakash, 
Vandenberg, & Buchak, 2015; Vandenberg, Prakash, & Buchak, 2014). 

The assumption of the model, that lateral variations in velocities, tem- 
peratures, and constituents are negligible, may be inappropriate for large 
waterbodies that exhibit significant lateral variations in water quality. 
Whether this assumption is met is often a judgment call by the user and 
depends in large part on the questions that are being addressed. 
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EFDC-EPA 

The Environmental Fluid Dynamics Code (EFDC) is a public domain, 
open source, surface water modelling system, which includes hydro- 
dynamic, sediment and contaminant, and eutrophication modules fully 
integrated in a single source code implementation. EFDC was originally 
developed at the Virginia Institute of Marine Science (VIMS) by Dr. John 
M. Hamrick (Hamrick, 1992) in 1988 and was later enhanced by Tetra 
Tech for the USEPA (Tetra Tech, 2007). The model is currently maintained 
by Tetra Tech with support from the USEPA. 

EFDC isa state-of-the-art hydrodynamic and water quality model that 
can be used to simulate aquatic systems in one, two, and three dimensions. 
EFDC uses stretched or sigma vertical coordinates and Cartesian or curvi- 
linear-orthogonal horizontal coordinates to represent the physical char- 
acteristics of a waterbody. For one-dimensional applications, an optional 
cross-section description can be used. Two horizontal grid generation and 
preprocessing tools, GEFDC (GridEFDC) and VOGG (Visual Orthogonal 
Grid Generator), are also available. Based on a semi-implicit conserv- 
ative finite volume solution scheme, EFDC’s hydrodynamic component 
solves three-dimensional, vertically hydrostatic, free surface, turbulent 
averaged equations of motion for a variable-density fluid. Dynamically 
coupled transport equations for turbulent kinetic energy, turbulent length 
scale, salinity, and temperature are also solved. The EFDC model allows 
for drying and wetting in shallow areas by a mass conservation scheme. 
Additional capabilities include simulation of shoreline movement by 
drying and wetting, hydraulic control structures, vegetation resistance, 
wave-current boundary layers, and wave-induced currents. 

The EFDC simulates multiple size classes of cohesive and non-co- 
hesive sediments. A sediment processes function library allows the user 
to choose from a wide range of currently accepted parameterizations for 
settling, deposition, resuspension, and bed load transport. The sediment 
bed is represented by multiple layers and includes several armoring rep- 
resentations for non-cohesive sediment and a finite strain consolidation 
formulation for dynamic simulation of bed layer thickness, void ratio, and 
pore water advection. The sediment transport component can operate in a 
morphological mode with full coupling, with the hydrodynamic compon- 
ent representing the dynamic evolution of bed topography. 
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Table 6. Summary of Selected Receiving Water Quality Models 


Model& Spatial Type of Simulated Pro- Simulated GUI Availability References 
Source Dimension Simulation cesses Parameters and Support 
CE-QUAL-W2 2D laterally Event and Hydrodynamic and Water surface, Yes Public domain, Berger & Wells, 
(Portland averaged limited water quality (all velocity, open source 2014; Cole & Wells, 
State Univ.) rivers, long-term pollutants except temperature, 2017; Shabani et 
estuaries, simulations toxics and metals) nutrients, al., 2017 
lakes, and for stratified and multiple algae, 
reservoirs; non-stratified zooplankton, 
includes systems periphyton, 
hydraulic macrophyte 
structures species, DO, 
PH, alkalinity, 
multiple 
CBOD, 
suspended 
solids, organic 
matters, and 
generic water 
quality groups 
EFDC (USEPA, 1D, 2D, and Event and Hydrodynamic, Temperature, Yes Public domain J. Craig et al., 
VIMS, Tetra 3D rivers, long-term sediment transport, cohesive and 2007; Hamrick, 
Tech, Inc.) estuaries, simulations water quality (eu- non-cohesive 1992: Hua & 
lakes, and trophication, sed- sediments, Zhang, 2017; Osmi, 
coastal waters iment diagenesis), COD/DO, Ishak, Kim, Azman, 
and hydraulic toxics (adsorption, nutrients, & Ramli, 2016; Seo, 


structures 


degradation) 


algae, salinity, 
metals, 

and other 
contaminants 


Sigdel, Kwon, & 
Lee, 2010; Tetra 
Tech, 2007 
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Table 6. (continued) 


EFDC-Plus 1D, 2D, and Event and Hydrodynamic, Temperature, Yes Proprietary DSI E. Cho, 

(DSI) 3D rivers, long-term sediment transport cohesive and for DSI multi- Arhonditsis, Khim, 
estuaries, simulations (added SEDZLJ non-cohesive thread version; Chung, & Heo, 
lakes, and approach), water sediments, public open 2016, 2017; Ji, 
coastal waters quality (eutrophi- COD/DO, for DSI single 2017; Shen et al., 
and hydraulic cation, sediment nutrients, thread version 2014 
structures diagenesis, rooted algae, salinity, 

plant and epiphyte), metals, 
toxics (adsorption, and other 
degradation, vola- contaminants 
tilization) 

HEC-RAS 2D overland, Event and 1D steady non- Bacteria, Yes Public domain Brunner, 2016; 

(USACE) 1D and 2D long-term uniform and temperature, Knebl, Yang, 
rivers, lakes, simulations unsteady flows sediments, Hutchison, & 
reservoirs, for WS profiles, BOD/DO, Maidment, 2005; 
and inundated floodways salinity Patel, Ramirez, 
floodplains, and floodplain Srivastava, Bray, 
and hydraulic determination, & Han, 2017; Wu & 
structures sediment and Fan, 2017; Xiong, 

limited water 2011 
quality. 2D 

capabilities are 

recently added 

MIKE 11 (DHI) 1D river Event and River hydraulics Sediment, Yes Proprietary DHI, 2017a; Liang 
reaches and long-term and sediment temperature, et al., 2015; 
hydraulic simulations transport; links to BOD/DO, Patro, Chatterjee, 
structures MIKE 21 for 1D and salinity, Mohanty, Singh, 

2D flood and MIKE nutrients & Raghuwanshi, 


SHE & ECOLAB 
for water quality 
simulations 


2009; Thompson, 
S@renson, Gavin, & 
Refsgaard, 2004; 
Zhao, Zhang, 
James, & Laing, 
2012 
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Table 6. (continued) 


Model& Spatial Type of Simulated Pro- Simulated GUI Availability References 
Source Dimension Simulation cesses Parameters and Support 
QUAL2K/ 1D rivers Event and Quasi-dynamic Sediment, Yes Public domain Brown & Barnwell, 
QUAL2Kw and streams long-term simulations with temperature, 1987; Chapra & 
(USEPA, divided into simulations steady (QUAL2k) BOD/DO, Pelletier, 2003; 
Washington sub-reaches or or non-steady salinity, pH Pelletier & Chapra, 
Department of computational state (QUAL2Kw) and alkalinity, 2005; Salvai & 
Ecology) elements state hydraulics, nutrients, Bezdan, 2008 

non-uniform steady algae 

flow, repeating ae 

diel conditions, periphyton, 

and water-quality pathogen 

kinetics. Has auto- 

calibration and 

uncertainty analysis 

capabilities 
WASP 1D,2D,and3D Event and With links to Bacteria, Yes Public domain Ambrose & Wool, 
(USEPA) rivers, lakes, long-term 1D, 2D, and 3D sediments, 2009; Ambrose, 

reservoirs, simulations hydrodynamic BOD/DO, Wool, Connolly, & 
estuaries, and models for dynamic nutrients, Schanz, 1988; Di 


coastal waters 


flow inputs, it 
simulates water 
temperature, three 
types of sediments, 
biochemical 
oxygen demand, 
sediment oxygen 
demand, dissolved 
oxygen, nitrogen, 
phosphorus, 
multiple species 

of algae, detritus, 
periphyton, organic 
toxicants, mercury 
and other metals, 
PH and alkalinity, 
and pathogens 


toxic organics, 
toxic metals, 
mercury, 
salinity, pH, 
and alkalinity 


Toro, Fitzpatrick, & 
Thomann, 1983; J. 
M. Johnston et al., 
2017; Z. Liu et al., 
2008 


The EFDC model includes a variable configuration eutrophication com- 
ponent for simulation of aquatic carbon, nitrogen, and phosphorus cycles. 
The kinetic processes included in the EFDC water quality model are de- 
rived from the US Army Corps of Engineers’ CE-QUAL-ICM water qual- 
ity model (Cerco & Cole, 1995), including sediment diagenesis. In contrast 
to earlier water quality models (such as WASP) (Ambrose & Wool, 2009), 
which use biochemical oxygen demand to represent oxygen demanding 
organic material, the EFDC water quality model is carbon-based and uses 
chemical oxygen demand (COD). The four algae species are represented in 
carbon units. The three organic carbon variables play an equivalent role 
to BOD. Organic carbon, nitrogen, and phosphorous can be represented 
by up to three reactive sub-classes, refractory particulate, and particulate 
and dissolved labile. In addition to the internal eutrophication model, 
EFDC can create hydrodynamic transport files formatted for WASP and 
CE-QUAL-ICM. 

EFDC can also represent the transport and fate of an arbitrary num- 
ber of contaminants, including metals and hydrophobic organics, sorbed 
to any of the sediment classes and dissolved and particulate organic car- 
bon using a three-phase equilibrium partitioning formulation. Dissolved 
and particulate organic carbon can be represented as independent state 
variables, and pollutants of concern can be fractionally assigned to any of 
the sediment classes. A contaminant processes function library allows the 
representation of various degradation and transformation processes. 

In addition to the grid generation tools, a windows-based model 
interface, EFDCView, which incorporates grid generation, pre-processing 
and post-processing tools, is available. The EFDC has been used for many 
modelling studies of rivers, lakes, estuaries, coastal regions, and wetlands 
internationally (Hua & Zhang, 2017; Huang, Falconer, & Lin, 2017; Osmi 
et al., 2016; Seo et al., 2010). Due to its range of applicability with respect to 
water body and pollutant types, EFDC has been a choice model for TMDL 
development in the United States, such as the Christina River (Merrill et 
al., 2002), Wissahickon Creek (Zou et al., 2006), Tenkiller Ferry Lake (P. 
M. Craig, 2006), Los Angeles Harbor (J. Craig et al., 2007), and Charles 
River (Peng et al., 2011). 
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EFDC-Plus (EFDC+) 
The ongoing evolution of the EFDC model has been application-driven by 
a diverse group of EFDC users in the academic, government, and private 
sectors. Since 2002, Dynamic Solutions International, LLC (DSI) has been 
steadily improving and enhancing the original EPA version of the EFDC 
code. The improvements to the EFDC code had become so extensive that 
in 2016 the DSI version of EFDC was renamed as the EFDCPlus or EFDC+ 
model (DSI, 2017; E. Cho et al., 2016; 2017; Ji, 2017; Shen et al., 2014). The 
EFDC+ and its associated powerful graphical user interface, EFDC_ 
Explorer, constitute the modelling system called the EFDC_Explorer 
Modelling System (EEMS), which supports integrated hydrodynamic, 
sediment, water quality, and toxics modelling for receiving water bodies. 
Compared with EFDC-EPA, EFDC+ has many key enhancements, for 
example: 


e Dynamic memory allocation: Dynamic memory 
allocation allows the user to use the same executable 
code for applications to different water bodies. 
Dynamic allocation eliminates the need to re-compile 
the EFDC code for different applications, because of 
different maximum array sizes required to specify the 
computational grid domain and time series input data 
sets. Dynamic allocation also helps prevent inadvertent 
errors and provides better traceability for source code 
development. 


e OpenMP - Multithreading: OpenMP provides vastly 
improved model run times. The Intel” OpenMP* 
Runtime Library binds OpenMP threads to physical 
processing units. Depending on the machine topology, 
application, and operating system, thread affinity can 
have a substantial impact on the application speed. 
EFDC+ typically produces run time up to four times 
faster on a six-core processor than the conventional 
single-threaded EFDC model. 


e Sigma-Zed layering: An improved version of the 
EFDC code helps deal with pressure gradient errors 
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that occur in simulations that have steep changes in 
bed elevation. The Sigma Zed code contrasts with 

the conventional EFDC code, which uses a sigma 
coordinate transformation in the vertical direction 

and uses the same number of layers for all cells in the 
domain. In the EFDC_SGZ model, the vertical layering 
scheme has been modified to allow for the number of 
layers to vary over the model domain. This approach 

is computationally efficient and provides significantly 
improved simulations of thermal stratification. 


Hydraulic structures: Equations governing hydraulic 
structures such as culverts, weirs, sluice gates, and 
orifices are implemented in EFDC+. This additional 
feature is different from the previous head lookup table 
used to describe the relationship between head and flow 
for a hydraulic structure. 


Enhanced heat exchange: Heat exchange options use 
equilibrium temperatures for the water and atmospheric 
interface and spatially-variable sediment bed 
temperatures. 


Ice formation and melt. 


Lagrangian particle tracking: This option is applicable to 
oil spill modelling and emergency response simulations. 


Improved/simplified external wave model linkage. 


SEDZL] toxics implementation: The SEDZLJ sedflume 
model developed by Sandia National Laboratories 
has been adopted and greatly enhanced in EFDC+ 

in addition to the original approach used in the EPA 
version of EFDC for sediment transport modelling. 
A whole new toxic modelling capability has been 
implemented. 


RPEM module: A Rooted Plant and Epiphyte Model 
(RPEM) has been incorporated into EFDC to better 
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simulate water quality interactions with submerged 
aquatic vegetation (epiphytic algae and macrophytes). 


e Internal wind wave generation: A wind-generated wave 
sub-model has been added to enable the computation 
of wind-generated wave bed shear stress on sediment 
resuspension and wave induced currents. 


e High frequency output: These are new output snapshot 
controls for targeting specific periods for high frequency 
output within the standard regular output frequency. 


e Restart/continuation run options: This application gives 
more options to run a model from the beginning or 
continue to run from the previous stop. 


e Streamlined code for quicker execution time. 


e Model linkages: This enhancement customizes linkage of 
model results for the Windows-based EFDC_Explorer 
graphical pre- and post-processor for EFDC+. 


e MHK linkage: Incorporating the Marine and Hydro- 
kinetic (MHK)-friendly module allows for the 
simulation of placement and potential effects of 
installing and operating turbines and wave energy 
converters in rivers, tidal channels, ocean currents, 
and other waterbodies. This code comes from Sandia 
National Laboratories modified Environmental Fluid 
Dynamics Code (SNL-EFDC) (Thanh, Grace & Carlton, 
2008). 


The aforementioned enhancements greatly improve the applicability of 
the EFDC model. 


HEC-RAS 

The Hydrologic Engineering Center River Analysis System (HEC-RAS) is 
a very powerful model for simulating one-dimensional, steady, non-uni- 
form hydraulics, and one-dimensional, two-dimensional, and combined 
one-/two-dimensional unsteady flow through a full network of open 


44 3.0 Modelling Approaches for Each Component 


channels, floodplains, and alluvial fans (Brunner, 2016). It evolved from 
the first FORTRAN version of HEC-2 released by the United States Army 
Corps of Engineers Hydrologic Engineering Center (HEC) in 1966. 

The HEC-RAS system contains several river analysis components for: 
(i) steady flow water surface profile computations, (ii) one- and two-di- 
mensional unsteady flow simulation, (iii) movable boundary sediment 
transport computations, and (iv) water quality analysis. In addition to 
these river analysis components, the system contains several hydraulic de- 
sign features that can be invoked once the basic water surface profiles are 
computed. The steady flow water surface profiles component is intended 
for calculating water surface profiles for steady, gradually varied flow. The 
system can handle a full network of channels, a dendritic system, or a 
single river reach. The unsteady flow component can be used to perform 
subcritical, supercritical, and mixed flow regime (subcritical, supercritic- 
al, hydraulic jumps, and drawdowns) calculations in the unsteady flow 
computations module. 

The sediment transport component of the modelling system is in- 
tended for the simulation of one-dimensional sediment transport/mov- 
able boundary calculations resulting from scour and deposition over 
moderate time periods (typically years, although applications to single 
flood events are possible). The sediment transport potential is comput- 
ed by grain size fraction, thereby allowing the simulation of hydraulic 
sorting and armoring. Major features include the ability to model a full 
network of streams, channel dredging, various levee and encroachment 
alternatives, and the use of several different equations for the computation 
of sediment transport. 

The water quality component of the modelling system is intended 
to allow the user to perform riverine water quality analyses. An advec- 
tion-dispersion module is included with this version of HEC-RAS, add- 
ing the capability to model water temperature. This new module uses 
the QUICKEST-ULTIMATE explicit numerical scheme to solve the 
one-dimensional advection-dispersion equation using a control volume 
approach with a fully implemented heat energy budget. The transport 
and fate of a limited set of water quality constituents is available in HEC- 
RAS, including: dissolved nitrogen (NO,-N, NO,-N, NH,-N, and Org-N), 
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dissolved phosphorus (PO,-P and Org-P), algae, dissolved oxygen (DO), 
and carbonaceous biological oxygen demand (CBOD). 

HEC-RAS includes a user-friendly interface and has a variety of data 
storage, management, and graphics and reporting components. A com- 
panion program, HEC-GeoRAS, provides tools and utilities for process- 
ing geospatial data in ArcGIS using a graphical user interface (GUI) for 
preparation of geometric data for import into HEC-RAS. 

The model is commonly used for simulating steady-flow water sur- 
face profiles or unsteady flow hydraulic and hydrodynamic simulations in 
support of hydraulic structure design, floodplain delineation, or floodway 
determination (Knebl et al., 2005; Patel et al., 2017; Xiong, 2011). Due to 
limited water quality capabilities, HEC-RAS is often used to generate hy- 
draulic field inputs for other water quality models (Hosseini et al., 2016; 
Wu & Fan, 2017). 


MIKE 11 

MIKE 11 is a one-dimensional fully dynamic model for simulating flows, 
sediment transport, and water quality in estuaries, rivers, irrigation chan- 
nels, and other water bodies (DHI, 2017a). The model is a part of the MIKE 
suite of water modelling software products developed by DHI Water and 
Environment. 

The hydrodynamic (HD) module is the nucleus of the model, which 
solves the vertically-integrated equations for conservation of continuity 
and momentum, i.e. the Saint Venant equations. Advanced computational 
modules are included for the description of flow over hydraulic structures, 
including possibilities to describe structure operation. The primary fea- 
ture of the MIKE 11 modelling system is the integrated modular structure 
with a variety of ad-on modules each simulating phenomenon related to 
river systems, including advection-dispersion, cohesive and non-cohesive 
sediment transport, and water quality. 

The advection-dispersion (AD) module is based on the one-dimen- 
sional equation of conservation of mass of dissolved or suspended ma- 
terial, i.e. the advection-dispersion equation. Non-cohesive standard and 
advanced cohesive sediment transport modules are part of the AD mod- 
ule. The module requires output from the hydrodynamic module, in time 
and space, as well as in terms of discharge and water level, cross-sectional 
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area, and hydraulic radius. The advection-dispersion equation is solved 
numerically using an implicit finite difference scheme which, in princi- 
ple, is unconditionally stable and has negligible numerical dispersion. In 
association with DHI ECOLAB water quality analysis and simulation of 
fate and transport in riverine systems, it can be used to develop TMDL for 
a variety of constituents (Cabrejo, 2011; Liang et al., 2015). MIKE 11 can 
be dynamically linked to DHI software, MIKE 21, to perform two-dimen- 
sional river and floodplain simulaions (Patro et al., 2009), and MIKE SHE 
for surface-groundwater interactions (Thompson et al., 2004; Zhao et al., 
2012). 

With its exceptional flexibility, speed, and user-friendly environment, 
MIKE 11 is an effective modelling tool to support detailed analysis, de- 
sign, management, and operation of channel systems. MIKE 11 has been 
used in numerous applications around the world for flood plain analysis 
and mapping, real-time flood, inflow and water quality forecasting, an- 
alysis and design of hydraulic structures, sediment transport, dredging 
impact and channel restoration alternative analysis, water quality an- 
alysis, issues related to TMDL and ecosystem restoration, and integrated 
groundwater and surface water analysis. In many aspects, MIKE 11 is very 
similar to HEC-RAS, with added benefits: direct linkage with the water- 
shed and groundwater flow components of MIKE SHE to allow integrated 
hydrologic, hydraulic, and hydrogeological modelling; a soft linkage to 
ECOLAB to allow water quality analysis and TMDL development; and 
a hard linkage to MIKE 21 (MIKE FLOOD) to allow a combination of 
one-dimensional and two-dimensional flood simulations. 

DHI introduced MIKE HYDRO River as the new-generation river 
modelling software and as a successor to MIKE 11 (DHI, 2017b). The re- 
lease of MIKE 2017 includes both MIKE 11 and MIKE HYDRO River. 
MIKE HYDRO River includes most features and add-ons available in 
MIKE 11. 


QUAL2K 

QUAL2K (or Q2K) (Chapra & Pelletier, 2003) is a one-dimensional, river 
and stream water quality model (Brown & Barnwell, 1987). The QUAL2K 
code is distributed by the USEPA. A variation of QUAL2K is distributed 
as QUAL2Kw (Pelletier & Chapra, 2005) by the Washington Department 
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of Ecology, and includes versions with auto-calibration and Monte Carlo 
simulation. 

The QUAL2K model is an extremely powerful model based on many 
of the same assumptions as QUAL2E (an earlier version of QUAL2K), such 
as a one-dimensional system with steady-state, non-uniform flows and 
hydraulics, which allows simulation of diel variations in water quality. It 
simulates a wide variety of conventional pollutants as well. Enhancements 
over QUAL2E include algorithms for slow and fast carbonaceous biochem- 
ical oxygen demand, periphyton and detritus (in addition to sediment 
diagenesis), pH and alkalinity, and other advanced features. The model 
input and output are in the form of user-friendly Excel spreadsheets, with 
underlying VBA routines to write and read files for use in a FORTRAN 
executable code. 

QUAL2K is applicable to waste-load allocation (WLA) and TMDL 
studies (Salvai & Bezdan, 2008) of rivers, streams, and some estuaries, 
using tidally averaged dispersion coefficients for conventional pollutants 
such as pathogens, nitrogen, phosphorus, dissolved oxygen, biochemical 
oxygen demand, sediment oxygen demand, phytoplankton, benthic algae, 
and pH. It is not applicable to toxics or metals and is limited to simulation 
of steady time-invariant flow, while time variations in water quality are 
only over diel cycles but constant otherwise. 

In the recent version 6 of QUAL2Kw, the model has been updated into 
a dynamic water quality model which simulates non-steady, non-uniform 
flow using kinematic wave flow routing. Continuous simulations can be 
run with time-varying boundary conditions for periods of up to one year, 
with the option to use repeating diel conditions with either steady or non- 
steady flows. 


Water Analysis Simulation Program (WASP) 

WASP is a general dynamic mass balance framework for contaminant fate 
and transport in surface water aquatic systems, including both the water 
column and the underlying benthos. The model has been continuously 
supported by the USEPA and enhanced since its original development in 
the 1980s (Ambrose & Wool, 2009; Ambrose et al., 1988; Di Toro et al., 
1983). 
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Based on the flexible compartment modelling approach, WASP can 
be applied in one, two, or three dimensions with advective and dispersive 
transport between discrete physical compartments or segments. WASP is 
designed to permit easy substitution of user-written routines into the pro- 
gram structure to form different water quality modules. 

The WASP code can simulate water temperature, three types of 
sediments, biochemical oxygen demand, sediment oxygen demand, dis- 
solved oxygen, nitrogen, phosphorus, multiple species of algae, detritus, 
periphyton, organic toxicants, metals, pH and alkalinity, and pathogens. 
The software includes a data preprocessor to format input datasets from 
simple ‘cut and paste’ to detailed queries from a database. A post-proces- 
sor (MOVEM) provides an efficient method for reviewing simulations 
with field data for calibration and confirmation testing. Simulations can 
be transferred to spreadsheets as *.CSV files and plotted or animated 
two-dimensionally. 

WASP is one of the most widely used water quality models in the 
United States and throughout the world. Because of the model’s capacity to 
handle multiple pollutant types, it has been widely applied in the develop- 
ment of TMDLs. WASP has the capability of linking with hydrodynamic 
and watershed models, which allows for multi-year analysis under varying 
meteorological and environmental conditions (J. M. Johnston et al., 2017; 
Z. Liu et al., 2008). 


3.2.2 Watershed Water Quality Models 
Most of the watershed models were either designed for hydrologic mod- 
elling to support soil and water management, or for general watershed 
and other resource management. Water quality modelling capabilities are 
limited or absent in most of the existing watershed models. Watershed 
models have been applied to predict non-point pollutant loadings through 
surface runoff and provide inputs to receiving water quality models. The 
review herein is simply a screening of existing watershed models capable 
of conducting daily and/or sub-daily water quality modelling for potential 
inclusion in integrated modelling system. 

A review of selected watershed water quality models is summarized 
in Table 7 using eight model characteristics/features, followed by brief de- 
scriptions about the source, capabilities, and applications of the models. 
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Watershed water quality models commonly require extensive input data, 
including topography, hydrography with reach geometry, LULC, soils, 
meteorology, agricultural practices (e.g., crop rotation, schedules, tillage 
practices, fertilizer applications, pesticide and herbicide applications). 


AnnAGNPS 

The Annualized Agricultural Non-Point Source Pollution Model 
(AnnAGNPS) is a watershed-scale continuous simulation model, which 
is an expansion of the capabilities developed in the single event mod- 
el AGNPS. The model is freely available from the US Department of 
Agriculture (USDA). 

On a daily timestep, AnnAGNPS simulates quantities of surface 
water, sediment, nutrients, and pesticides which are leaving the land areas, 
as well as their subsequent travel through the watershed. In the model, a 
watershed is subdivided into homogenous land areas with respect to soil 
type, land use, land management, and climate. Areas can be of any shape 
including hydrologically-based or square grid. The soil profile is divided 
into two layers. The top 200 mm is used as a tillage layer whose properties 
can change (bulk density, etc.). The remaining soil profile comprises the 
second layer whose properties remain static. A daily soil moisture water 
budget includes applied water (rainfall, irrigation, and snowmelt), runoff, 
evapotranspiration, and percolation. Runoff is calculated using the Soil 
Conservation Service runoff curve number equation but is modified if a 
frozen, shallow surface soil layer exists. Curve numbers are modified daily 
based upon tillage operations, soil moisture, and crop stage. 

Overland erosion of sediment is determined using the revised 
Universal Soil Loss Equation (RUSLE). Special components are included 
to handle concentrated sources of nutrients (feedlots and point sources), 
ephemeral gully sources, concentrated sediment sources (classical gullies), 
added water (irrigation), and the impacts of riparian buffers and wetlands. 
The model partitions soluble nutrients and pesticides between surface 
runoff and infiltration. Soluble nutrients from feedlots are also transport- 
ed with runoff. Sediment-transported nutrients and pesticides are also 
determined. The sediment generated from the land areas and gullies is 
subdivided into particle size classes (clay, silt, sand, small aggregate, and 
large aggregate) before being added to the stream system. Particle sizes 
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Table 7. Summary of Selected Watershed Water Quality Models 


Model & Spatial Rep- Type of Simulated Simulated Management BMP GUI Availability References 
Source resentation Simulation Processes Parameters Operations Measures and Support 
AnnAGNP Any Event and Hydrology, Surface Irrigation, Agri- Yes Public Bingner, 
(USDA- shapes of continuous snowmelt, runoff, pumping cultural domain Theurer, & 
ARS) homogenous (daily) plant subsurface BMPs Yuan, 2015; 
land areas, growth, land lateral flow, Yasarer et 
channels management, sediments, al., 2017 
erosion, nutrients, 
pollutant chemical 
loadings, fate, oxygen de- 
and transport mand, and 
pesticides 
HEC-HMS Sub- Event and Precipita- Surface and Reservoir No Cin- Yes Public Pak, 
(US Army watersheds, continuous tion, snow subsur- operations directly domain Fleming, 
Corps of reaches, and accumulation face flows, reflected Scharffen- 
Engineers) junctions and melting, snow melt, vis the berg, 
direct runoff sediment, topo- Gibson, 
(overland nitrogen, graphic & Brau- 
flow and phospho- factor er, 2015; 
interflow), rus, algae, describ- Scharffen- 
baseflow, dissolved ing the berg, 2016 
flow routing, oxygen influence 
infiltration, (DO), and of plant 
evapotranspi- carbo- cover on 
ration naceous surface 
biological erosion) 
oxygen 
demand 
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Table 7. ( continued) 


Model & Spatial Rep- Type of Simulated Simulated Management BMP GUI Availability References 
Source resentation Simulation Processes Parameters Operations Measures and Support 
HSPF Catchments Event and Hydrology, Surface and  Land-use Yes, Yes Public ASCE, 
(USEPA with continuous snowmelt, subsurface management moderate domain 2017; 
and USGS) pervious and erosion, flows, snow practices, level of Bicknell, 
impervious pollutant melt, con- and water analysis; Imhoff, 
areas, loadings, fate, servatives, management some lim- Kittle, 
channels, and and transport sediment, operations itations. Jobes, & 
reservoirs tempera- Explicit Donigian, 
ture, DO, BMP 2005 
biochemical represen- 
oxygen de- tation in 
mand, pH, version 
ammonia, 12.4 to 
nitrite-ni- release 
trate, organ- 
ic nitrogen, 
orthophos- 
phate, 
organic 
phospho- 
rous, phy- 
toplankton, 
zooplank- 
ton, fecal 
coliforms, 
and pesti- 
cides 


Table 7. ( continued) 


LSPC Catchments 

(Tetra Tech with pervious 

Inc.) & impervious 
areas, 
channels, and 
reservoirs 
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Hydrology, 
snowmelt, 
erosion, 
pollutant 
loadings, fate, 
and transport 


Surface and 
subsurface 
flows, snow 
melt, con- 
servatives, 
sediment, 
tempera- 
ture, DO, 
biochemical 
oxygen de- 
mand, pH, 
ammonia, 
nitrite-ni- 
trate, organ- 
ic nitrogen, 
orthophos- 
phate, 
organic 
phospho- 
rous, phy- 
toplankton, 
zooplank- 
ton, fecal 
coliforms, 
and pesti- 
cides 


Land-use 
management 
practices, 
and water 
management 
operations 


Yes, Yes 
moderate 

level of 
analysis; 

some lim- 
itations 


Public 
domain 


Tetra Tech, 
2009 


€S 
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Table 7. ( continued) 


Model & Spatial Rep- Type of Simulated Simulated Management BMP GUI Availability References 
Source resentation Simulation Processes Parameters Operations Measures and Support 
MIKE SHE 2D overland, Event and Precipitation, Surface and Irrigation, No, Yes Proprietary DHI, 2017c, 
(DHI 1D channels, continuous snow subsurface plumping, needs 2017d; 
1D-unsatu- accumulation flows, water control add-on Graham 
rated, and and melting, sediments, structures modules & Butts, 
3D-saturated overland and nutrients, 2005 
zones channel flow, pesticides, 
unsaturated and other 
zone, water 
saturated quality 
zone, parameters 
exchanges with link to 
between ECOLAB 


aquifers and 
rivers, crop 
growth and 
nitrogen 
processes in 
root zone, 
geochemical 
processes, 


advection and 


dispersion of 
solutes, soil 
erosion 
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Table 7. ( continued) 


SWAT Subwa- Event and Precipitation, | Weather, Irrigation, Agri- Yes Public Arnold, 

(USDA- tersheds, continu- snowmelt, surface and plumping, cultural domain Williams, 

ARS) channels, and ous (daily, surface run- subsurface reservoir BMPs Srinivasan, 

ponds sub-daily) off, interflow, flows, sed- King, & 

ground- iment, soil Griggs, 
water flow, tempera- 1994; 
infiltration, ture, crop Neitsch, 
percolation, growth, Arnold, 
evapotrans- nutrients, Kiniry, & 
piration, soil pesticides, Williams, 
temperature and agricul- 2011; 

i tural man- Pignotti, 
crop growth, agements Rathjens, 
soil erosion, Cibin, 
nutrients, Chaubey, & 
pesticides Crawford, 

2017; Sood 
& Ritter, 
2010; 
White & 
King, 2003 
SWMM Catchments, Event and Surface and Surface Stormwater Low Yes Public Alamdari, 
(USEPA) stream continuous subsurface runoff, and sewer impact domain Sample, 
& sewer flows, urban subsurface water develop- (proprietary Steinberg, 
network, & storm and flow, management ment PCSWMM, Ross, & 
ponds sanitary dynamic XPSWMM) Easton, 
sewer flows, flow routing 2017; 
sediment, in stream, Bhowmick, 
water quality subsurface Irvine, & 
drainage Jindal, 
network 2017; Ricks, 
and sanitary 2015; 
sanitary Rossman, 
sewers, and 2015 


loadings for 
up to ten 
pollutants 


are routed separately in the stream reaches. A Windows-based interface 
provides capabilities to subdivide the watershed into hydrologically de- 
rived cells, an input editor assisting in preparation of AnnAGNPS input 
data, and a processor that can calculate output loads at any point in the 
watershed. 

AnnAGNPS can be used to evaluate non-point source pollution from 
agricultural watersheds and to compare the effects of implementing vari- 
ous best practices over time within the watershed (Yasarer et al., 2017). 
Cropping and tillage systems for sheet and rill and ephemeral gully ero- 
sion, fertilizer, pesticide, irrigation application rates, point source loads, 
feedlot management, riparian buffers, and wetland management can be 
evaluated. However, there are several limitations that might restrict the 
application of the model for certain purposs: (i) all runoff and associated 
sediment, nutrient, and pesticide loads for a single daily event are rout- 
ed to the watershed outlet before the next day simulation; (ii) there is no 
tracking of nutrients and pesticides attached to sediment deposited in 
stream reaches from one event to the next event; and (iii) point sources 
are limited to constant loading rates (water and nutrients) for the entire 
simulation period. 


HEC-HMS 

The Hydrologic Modelling System (HEC-HMS) is designed by the United 
States Army Corps of Engineers Hydrologic Engineering Center to simu- 
late the complete hydrologic processes of dendritic watershed systems 
(Scharffenberg, 2016). The model is available to the public, with versions 
on Windows, Linux, and Solaris platforms. 

The hydrologic forecasts are based on a physical description of water- 
sheds, obtainable from geographic information systems combined with 
meteorological information and model parameters. The physical rep- 
resentation of a watershed is accomplished with a basin model formulated 
as a dendritic network which connects hydrologic elements such as sub- 
basin, reach, junction, reservoir, diversion, source, and sink. Computation 
for runoff processes proceeds from upstream elements in a downstream 
direction. A variety of methods are available for simulating hydrological 
processes, such as: (i) infiltration losses — initial constant, Soil Conservation 
Service (SCS) curve number, exponential, Green Ampt, Smith Parlange, 
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deficient and constant, and soil moisture accounting methods; (ii) sur- 
face runoff - Clark, modified Clark, Snyder, SCS and user-specified unit 
hydrograph methods, and the kinematic wave method; (iii) baseflow - re- 
cession, bounded recession, constant monthly, linear reservoir, and non- 
linear Boussinesq methods; (iv) open channel flow - kinematic wave, lag, 
modified plus, Muskingum, Muskingum-Cunge, and Straddle Stagger 
methods; (v) water impoundments - a user-entered storage-discharge re- 
lationship, pumps, and physical spillway and outlet structures; and (vi) 
diversion structures — user-specified function, lateral weir, pump station, 
and observed diversion flows. 

Sediment and water quality are simulated in an optional component 
in the basin model of the HEC-HMS. Sediment yields are estimated from 
land surface erosion and channel and reservoir transport. Sediment trans- 
port simulations define the non-cohesive sediment carrying capacity of 
flow in channels and reservoirs by grain-size distribution. The basin mod- 
el solves the advection-diffusion equation using the QUICKEST scheme. 
The channel reaches and the reservoir’s mass balances simulate nitrogen 
(as organic, ammonia, nitrite, and nitrate), phosphorus, algae, dissolved 
oxygen (DO), and carbonaceous biological oxygen demand. 

While the HEC-HMS model has been applied for hydrologic and flood 
studies, only few applications for water quality purposes are reported in 
the literature. The HEC-HMS model has been applied to simulate soil ero- 
sion and sediment transport in the Upper North Bosque River Watershed 
(UNBRW).The HEC-HMS results matched the observed TSS at five gauge 
locations across the UNBRW (<1% error at all gauges) during model cali- 
bration, and maintained modest residuals (-31 to 12% error) during the 
validation period (Pak et al., 2015). 


HSPF 

The Hydrologic Simulation Program-Fortran (HSPF) model developed 
in the 1960s (supported by both the USEPA and the USGS) has a long 
history of development and applications. In the 1990s, HSPF was selected 
as the core watershed model in the BASINS (Better Assessment Science 
Integrating Point and Non-point Sources) modelling system (K. Borah 
& Bera, 2004). USEPA is expected to release version 12.4 with a number 
of enhancements including explicit Best Management Practice (BMP) 
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representation and modelling, dynamic wave channel flow routing, and 
wetland modelling capabilities (ASCE 2017). 

HSPF is a process-based semi-distributed watershed model capable 
of simulating a single event, as well as continuous water quantity and 
quality, in urban and rural watersheds. It uses sub-basins as hydrologic 
response units and represents the landscape as pervious, impervious, and 
water body reach segments (Bicknell et al., 2005). The model has three 
basic modules: (i) PERLND (Pervious Land Upland Loading Module) 
represents hydrologic and water quality processes that are specific to 
pervious surfaces, (ii) IMPLND (Impervious Land Module) represents 
processes specific to impervious surfaces, and (iii) RCHRES (Reservoir 
Routing Module), which is a one-dimensional stream model that serves 
as the receiving water model. The PERLND module has upper, lower, 
and active groundwater zones. It uses simple storage-based equations for 
one-dimensional flow routing. The HSPF uses meteorological input time 
series data and computes hydrology and water quality time series data. 
The model simulates interception, soil moisture, surface runoff, interflow, 
base flow, snowpack depth and water content, snowmelt, evapotranspir- 
ation, groundwater recharge (flux to deep aquifer), dissolved oxygen, 
biochemical oxygen demand, temperature, pesticides, conservative con- 
stituents, fecal coliform, sediment detachment and transport, sediment 
routing by particle size, channel routing, reservoir routing, constituent 
routing, pH, total nitrogen, total phosphorus, ammonia, nitrite-nitrate, 
organic nitrogen, orthophosphate, organic phosphorus, phytoplankton, 
and zooplankton. The model operates on any timestep from one minute 
to one day. Because inputs for many parameters are needed to character- 
ize the lumped hydrologic response units at sub-basin level, calibrating a 
comprehensive HSPF model is generally a challenging task. 

HSPF has been used for a variety of applications such as assessing 
the impact of land-use change, reservoir operations, point or non-point 
source pollution management, flow diversion and water withdrawal im- 
pact assessment studies, and setting TMDLs for water quality-impaired 
water bodies (ASCE, 2017). It can model areas from one square meter to 
thousands of square kilometers at user-specified timesteps. These capabil- 
ities enable model users to evaluate the impact of BMPs over many years 
through wet and dry cycles and the dynamic evaluation of pollutant 
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loadings to receiving water bodies. Although HSPF is a comprehensive 
and highly flexible model, the current version of the model has limited 
abilities to explicitly represent and simulate BMPs such as detention basins 
and infiltration trenches. 


LSPC 

The Loading Simulation Program in C++ (LSPC) is a watershed modeling 
system that includes streamlined HSPF algorithms for simulating hydrol- 
ogy, sediment, and general water quality on land, as well as a simplified 
stream fate and transport model (Tetra Tech, 2009). LSPC was originally 
developed by Tetra Tech for the USEPA Region 3 in the early 2000s to 
determine watershed-scale TMDLs. The code is publicly available from 
the USEPA. 

The primary difference between the LSPC and the HSPF is the pro- 
gramming architecture. The LSPC uses C++ to use common data manage- 
ment software (i.e., Microsoft Access) and to avoid inherent limits on data 
array size and spatial and temporal resolution. To streamline interpreta- 
tion for each simulation, the LSPC automatically generates comprehensive 
sub-watershed files for all land-layers, reaches, and simulated modules 
using hourly or daily intervals. The software also calculates the TMDL 
and allocates the source reductions. 

Because of its C++ programming architecture, which has no inherent 
limits on array size and spatial/temporal resolution associated with model 
setup, LSPC overcomes many of the difficulties experienced with large- 
scale watershed simulation. LSPC is frequently used for watershed appli- 
cations and can be readily linked to receiving water models such as EFDC, 
WASP, and CE-QUAL-W2 for complex waterbody TMDL development. 


MIKE SHE 
MIKE SHE is a distributed and physically based integrated surface and 
subsurface hydrological and water quality modelling system for simulat- 
ing the entire land phase of a hydrologic cycle. It is a proprietary model 
developed and maintained by the Danish Hydraulic Institute (DHI, 2017c, 
2017d). 

In the model, the study area is divided into polygons based on land 
use, soil type, and precipitation region. The polygons are then assigned 
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identification numbers. Model input files can be generated by overlaying 
the model input parameters with a grid network. Most of the data prepar- 
ation and model set-up can be completed using an external GIS software 
or MIKE SHE'’s built-in graphic preprocessor. The system has no limita- 
tions regarding watershed size. It can be implemented as a single event or 
as a long-term continuous simulation model using different time inter- 
vals for processes in different zones (overland hydrology, river hydraulics, 
groundwater, and water quality). 

MIKE SHE simulates the entire land phase of the hydrologic cycle 
from rainfall to stream flow and various flow processes, such as evapo- 
transpiration from vegetated land cover and evaporation from water bod- 
ies, overland flow, infiltration into soils, unsaturated zone flow, ground- 
water flow, and interaction between groundwater and surface water bodies 
(interflow and base flow). MIKE SHE offers several different approaches 
for hydrologic processes ranging from simple, lumped, and conceptual 
approaches, to advanced, distributed, and physically-based approaches. 
With distributed physically based approaches, overland flow is routed 
using a two-dimensional finite difference solution of the diffusive wave 
approximation of the Saint Venant equation. Flow in the unsaturated zone 
is solved using a one-dimensional finite difference solution of the Richards 
equation, and flow in the saturated zone is solved using a three-dimen- 
sional finite difference solution for Darcy’s law. MIKE SHE can be coupled 
with MIKE 11, a one-dimensional hydraulic model based on the one-di- 
mensional solution for the Saint Venant equations, to conduct integrated 
watershed and receiving water modelling. The model also simulates water 
use and management operations including irrigation systems, pumping 
wells, and various water control structures. A variety of agricultural prac- 
tices and environmental protection alternatives may be evaluated using 
other add-on modules, such as MIKE SHE DAISY, the crop yield and 
nitrogen consumption module (Booty & Benoy, 2009). 

A generic ecological modelling tool called MIKE ECO Lab can be 
coupled to a MIKE SHE hydrologic model to allow for the representation 
of a range of water quality and ecological processes with respect to the 
river, surface water, soil, and groundwater systems. MIKE ECO Lab relies 
on other models to calculate flow and transport processes. With existing 
or customized MIKE ECO Lab water quality templates, MIKE SHE can be 
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applied to simulate the fate and transport of different substances (such as 
sediments, nutrients, and pesticides) across all hydrologic and hydraulic 
model components. MIKE SHE can simulate fully integrated solute trans- 
port between surface water and the subsurface, including decay, sorption, 
precipitation, and selective plant uptake. More complex, multispecies, and 
kinetic reactions comprising all aspects of eco-hydrology can also be set 
up with MIKE ECO Lab. 

As one of the first commercially available, distributed, and physical- 
ly-based hydrologic models (El-Nasr et al., 2005), MIKE SHE has been 
used in a wide range of research and application projects, including surface 
water and groundwater quality assessment and remediation, impacts of 
land use and climate changes on long term water availability and quality, 
and impacts of agricultural management practices (e.g., irrigation, drain- 
age, sediments, nutrients, and pesticides) (Graham & Butts, 2005). In the 
United States, most of the applications have been in Florida, where there 
are strong interactions between surface water and groundwater aquifers 
because of high water table conditions and extensive lakes and wetlands. 
For example, MIKE SHE and MIKE ECO Lab have been used to establish 
a water quality model for TMDL development in the city of Kissimmee in 
Central Florida. The model was calibrated with historical data and then 
used to identify significant pollutant load areas for the implementation of 
specific corrective measures to reduce quantities of loadings entering the 
receiving water bodies. 


SWAT 
The Soil and Water Assessment Tool (SWAT) is a semi-distributed, con- 
ceptual and continuous-time river basin or watershed-scale model de- 
veloped at the United States Department of Agriculture Agricultural 
Research Service (Arnold et al., 1994; Neitsch et al., 2011). The model is 
open source, has multiple geographic information system interfaces such 
as ArcSWAT and QSWAT for creating input files, and has user-friendly 
tools for model calibration. SWAT is included in the USEPA BASINS for 
non-point source simulations on agricultural lands. 

SWAT was developed to predict the impacts of land management prac- 
tices on water, sediment, and agricultural chemical yields in large complex 
watersheds with varying soils, land use, and management conditions over 
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long periods of time. In SWAT, a watershed is portioned into sub-water- 
sheds, and each sub-watershed is further divided into a number of hydro- 
logic response units (HRU). An HRU represents a lumped land area 
having unique soil, land cover, and land management characteristics. A 
channel network can be delineated at a chosen scale of interest with a 
watershed containing at least one main channel or reach and a tributary 
channel. Ponds, wetlands or reservoirs, and point sources can be added as 
additional subunits. 

SWAT is capable of simulating physical, chemical, biological, and 
physiological processes of watersheds. The model simulates water flow 
and water quality in uplands, in subsurface water, in stream channels, and 
in open water bodies (ponds, wetlands, and impoundments). SWAT is a 
process-based model which conserves water and constituent mass. For 
example, the land phase of the hydrologic cycle is based on the water bal- 
ance equation. The model uses a modification of the SCS curve number 
method or the Green and Ampt infiltration equation to compute surface 
runoff volume for each hydrologic response unit. Evapotranspiration is 
calculated from potential evapotranspiration, which is estimated by the 
Penman-Monteith, Priestly-Taylor, or Hargreaves method. Peak runoff 
is estimated using a modification of the rational method. Flow is routed 
through channels using either a variable storage coefficient method or the 
Muskingum routing method. The water balance for the simulated shal- 
low aquifer is expressed in terms of a linear flow balance equation, which 
is solved to yield a simple algebraic relationship for base flow. Similarly, 
mass conservation of snowmelt, sediment, nutrients, carbon, bacteria, 
pesticides, and dissolved oxygen are all formulated in terms of simple al- 
gebraic relationships akin to the explicit finite difference approximation of 
ordinary differential equations. A mass balance in vegetative filter strips, 
grassed waterways, wetlands, ponds, and impoundments or reservoirs is 
described by a similar numerical approach. As a management tool, SWAT 
simulates crop management, conservation, and agricultural best manage- 
ment practices, as well as water management. It can also be adapted to 
simulate management practices that are not explicitly represented in the 
model. 

Due to its open source nature and an active development commun- 
ity, SWAT is constantly being improved and augmented with new process 
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representations since it was created in the early 1990s. The model is best 
suited for long-term continuous applications. Most of the applications of 
SWAT have been on a daily timestep, although a recent addition to the 
model includes the Green and Ampt infiltration equation, using rainfall 
input at any time increment and channel routing at an hourly timestep. 
SWAT is ideally suited for addressing a wide array of issues related to 
climate change, land use change, bioenergy crops, blue and green water 
availability, sediment transport, nutrient cycle and contaminant loads, 
BMPs, and TMDL studies (Sood & Ritter, 2010; White & King, 2003). 
While the HRU approach provides a simple, computationally effi- 
cient framework, processes modelled on HRUs are lumped and therefore 
spatially disconnected, as they are routed directly to sub-basin outlets. 
This was identified as a key weakness of the model (Douglas-Mankin, 
Srinivasan, & Arnold, 2010). This lack of definition of landscape position 
makes implementation of spatially targeted management measures diffi- 
cult to incorporate into the model. To overcome the spatial limitations of 
the HRU approach, a grid-based version of the SWAT model, SWATgrid 
(Rathjens et al., 2015), was developed to perform landscape simulations on 
a regularized grid by employing a modified landscape routing algorithm. 
However, SWATgrid remains largely untested, with little understanding of 
the impact of user-defined model spatial resolution (Pignotti et al., 2017). 


SWMM 

The Storm Water Management Model (SWMM) is a comprehensive 
hydrologic and hydraulic model used for single event or long-term (con- 
tinuous) simulation of runoff quantity and quality primarily from urban 
areas. SWMM was originally developed by the USEPA in 1971 and has 
undergone several major upgrades (Rossman, 2015). SWMM is includ- 
ed in the BASINS and is publicly available. Two widely used proprietary 
software packages derived from SWMM are PCSWMM (Bhowmick et al., 
2017) and XP-SWMM (Ricks, 2015), which have additional graphic user 
interface functions and capabilities. 

The EPA SWMM has been widely used throughout the world for plan- 
ning, analysis, and design related to storm water runoff, combined sewers, 
sanitary sewers, and other drainage systems in urban areas, with many 
applications in non-urban areas as well. The current edition, version 5.1, 
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is a complete rewrite of the previous release and provides an integrated 
environment for editing data, running hydrologic, hydraulic, and water 
quality simulations. The runoff component of SWMM operates on a col- 
lection of sub-catchment areas that receive precipitation and generate 
runoff and pollutant loads. The routing portion of SWMM transports this 
runoff through a system of pipes, channels, storage/treatment devices, 
pumps, and regulators. SWMM accounts for various hydrologic process- 
es that produce runoff from urban areas, including nonlinear reservoir 
routing of overland flow, and capture and retention of rainfall/runoff with 
various types of low impact development practices. SWMM also contains 
a flexible set of hydraulic modelling capabilities used to route runoff and 
external inflows through a drainage system network of pipes, channels, 
storage/treatment units and diversion structures. SWMM can also esti- 
mate pollutant loads associated with runoff. The following processes can 
be modelled for any number of user-defined water quality constituents 


e dry weather pollutant build-up over different land uses, 


e pollutant wash-off from specific land uses during storm 
events, 


e direct contribution of rainfall deposition, 
e reduction in dry-weather build-up due to street cleaning, 
e reduction in wash-off load due to BMPs, 


e entry of dry weather sanitary flows and user-specified 
external inflows at any point in the drainage system, 


e routing of water quality constituents through the 
drainage system, and 


e reduction in constituent concentration through 
treatment in storage units or by natural processes in 
pipes and channels. 


The SWMM software has been used in thousands of studies worldwide, 
including storm water master planning, sewer master planning, flood- 
plain management, water quality and TMDL modelling, and BMP and 
LID evaluations (Alamdari et al., 2017). 
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3.3 Groundwater models 


Besides groundwater quantity, the quality of groundwater is equally im- 
portant as a wide variety of contaminants can be found in groundwater, 
both organic and inorganic, including synthetic organic chemicals, 
hydrocarbons, inorganic cations and anions, pathogens, and radionu- 
clides. The groundwater models focus on changes in storage and fluxes 
within the saturated zone and can be classified into physical, analogue, 
and mathematical models to simulate groundwater movement and con- 
taminant transport. Physical models can be developed in the laboratory 
to study specific problems of groundwater flow or contaminant transport. 
Analogue models are based on equations (e.g., Ohm’s law/Darcy’s law) 
which describes groundwater flow in isotropic homogenous porous media. 
Mathematical models rely on groundwater flow (differential) equations 
which can often be solved only by approximate methods using a numer- 
ical method. The most widely used numerical methods are finite element 
and finite difference methods. With the advancements in computing tech- 
nology, sophisticated groundwater models have been developed that can 
be interfaced with GIS or coupled with other models. A selective list of 
popular groundwater numerical models is presented in Table 8. 


MODFLOW 

MODFLOW (Modular Three-Dimensional Finite-Difference Groundwater 
Flow) was developed by the United States Geological Survey (USGS) 
for simulating steady-state and/or transient saturated groundwater flow 
under confined and unconfined conditions (Harbaugh, 2005). The model 
and source code are publicly available on the USGS website, along with 
many customized versions and utility programs. 

The model simulates three-dimensional groundwater flow through a 
porous medium, using a finite difference method originally documented 
by McDonald and Harbaugh (McDonald & Harbaugh, 1984). The model 
domain is discretized into rectangular grid cells, and spatial derivatives are 
approximated based on the difference between the values of groundwater 
head at neighboring nodes and the spatial distance between the nodes. 
The groundwater flow can then be reconstructed based on the potentio- 
metric heads. Numerous MODFLOW versions have been developed as a 
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result of growing interest in surface and groundwater interactions, solute 
transport, and saltwater intrusion (Christian D. Langevin et al., 2017). The 
customized packages include enhanced capabilities to simulate processes 
related to evapotranspiration, rivers, lakes, and multi-node wells (Jones & 
Mendoza, 2012). Public domain and commercial Graphical user interfaces 
(GUID are available for MODFLOW set up, execution and results post-pro- 
cessing, including the MODELMUSE maintained by the USGS. 

MODFLOW is the most widely used modelling tool in the world for 
simulating groundwater flow, with numerous applications by studies on 
surface and groundwater interactions (Barthel & Banzhaf, 2016; Golden et 
al., 2014; Guzman et al., 2015), climate change impacts (Chunn, Faramarzi, 
Smerdon, & Alessi, 2019), solute transport (Zhang et al., 2013),including 
Alberta Oil Sands environmental impacts assessment ((R. Thompson, 
Mooder, Conlan, & Cheema, 2011). 


MT3D-USGS 

MT3D-USGS (Bedekar et al., 2016) is a USGS updated release of the 
groundwater solute transport code MT3DMS (Modular Transport, 
3-Dimensional, Multi-Species model) version 5.3 (Zheng & Wang, 1999). 
MT3D-USGS includes new transport modelling capabilities to accommo- 
date flow terms calculated by MODFLOW packages that were previously 
unsupported by MT3DMS, and to provide greater flexibility in simulation 
of solute transport and reactive solute transport. MT3D-USGS is available 
in the public domain. 

MT3D-USGS uses simulated hydraulic heads, intercell flows, and 
source and (or) sink terms from the MODFLOW output in the solution of 
the advection dispersion equation. MT3D-USGS capabilities and features 
include 


e unsaturated-zone transport, 


e transport within streams and lakes, including solute 
exchange with connected groundwater, 


e capability to route solute through dry cells that 
may occur in the Newton-Raphson formulation of 
MODFLOW (that is, MODFLOW-NWT), 
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e chemical reaction option that includes the ability to 
simulate interspecies reactions and parent-daughter 
chain reactions, 


e pump-and-treat recirculation that enables the simulation 
of dynamic recirculation with or without treatment 
for combinations of wells that are represented in the 
flow model, mimicking the above-ground treatment of 
extracted water, 


e reformulation of the treatment of transient mass storage 
to improve conservation of mass and yield solutions for 
better agreement with analytical benchmarks, 


e separate specification of partitioning coefficient (Kd) for 
mobile and immobile domains, 


e capability to assign prescribed concentrations to the top- 
most active layer, 


e ability to ignore cross-dispersion terms, and 


e ability to specify an absolute minimum thickness rather 
than the default percentage minimum thickness in dry- 
cell circumstances. 


MT3DMS has been an industry standard, accepted by practitioners and 
researchers and applied in thousands of studies worldwide (Ghoraba, 
Zyedan, & Rashwan, 2013; H. Zhang, Xu, & Hiscock, 2013). Like MT3DMS, 
MT3D-USGS is designed as a generalized groundwater solute transport 
code for use with any block-centered finite-difference groundwater flow 
model such as MODFLOW. MT3D-USGS can be used to simulate changes 
in concentrations of contaminants in groundwater considering advection, 
dispersion, and chemical reactions. The chemical reaction package options 
available in the model include equilibrium-controlled linear or nonlinear 
sorption, first-order irreversible decay or biodegradation, interspecies re- 
actions, and parent-daughter chain reactions. 
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SUTRA 

SUTRA (Saturated-Unsaturated Transport) is a finite element simulation 
model for two-dimensional or three-dimensional saturated-unsatur- 
ated, fluid-density-dependent groundwater flow with energy transport 
or chemically-reactive single-species solute transport model (C. I. Voss & 
Provost, 2010). The original version of SUTRA was released in 1984, and 
the latest version, SUTRA 3.0, was released in 2019. 

The model employs a two-dimensional or three-dimensional fi- 
nite-element and finite-difference method to approximate the governing 
equations that describe the two interdependent processes that are simu- 
lated: fluid-density-dependent saturated or unsaturated groundwater flow 
and either (a) transport of a solute in the groundwater, in which the solute 
may be subject to equilibrium adsorption on the porous matrix and both 
first-order and zero-order production or decay, or (b) transport of thermal 
energy in the groundwater and solid matrix of the aquifer. SUTRA tracks 
the transport of either solute mass or energy in flowing groundwater 
through a unified equation, which represents the transport of either solute 
or energy. Solute transport is simulated through a numerical solution of 
a solute mass balance equation where the solute concentration may affect 
fluid density. The single solute species may be transported conservatively, 
or it may undergo equilibrium sorption (through linear, Freundlich, or 
Langmuir isotherms). In addition, the solute may be produced or decayed 
through first-order or zero-order processes. Energy transport is simulated 
through a numerical solution of an energy balance equation. The solid 
grains of the aquifer matrix and fluid are locally assumed to have equal 
temperature, which may affect fluid density and viscosity. 

As the primary calculated result, SUTRA provides fluid pressures 
and either solute concentrations or temperatures, as they vary with time, 
everywhere in the simulated subsurface system. SutraGUI is a public 
domain computer program designed to run with the proprietary Argus 
ONE package, which provides two-dimensional Geographic Information 
System (GIS) and meshing support (Winston & Voss, 2004). 

SUTRA’s modular design allows straightforward modifications to the 
code. Eventual modifications, for example, are the addition of non-equi- 
librium sorption (such as two-site models), equilibrium chemical reactions 
or chemical kinetics, or the addition of overburden and underburden heat 
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loss functions, a wellbore model, or confining bed leakage. The USGS’s 
SUTRA code is the most widely used simulator for seawater intrusion and 
other variable density groundwater flow problems based on solute trans- 
port or heat transport (C. Voss, 1999). It has also been widely used for 
many other types of problems (Tsanis, 2006). 


HGS 
HydroGeoSphere (HGS) is a three-dimensional control-volume finite-ele- 
ment simulator which is designed to simulate the entire terrestrial por- 
tion of the hydrologic cycle based on a rigorous conceptualization of the 
hydrologic system consisting of surface and subsurface flow regimes in 
fractured or unfractured porous media (Aquanty, 2015). Originally, it was 
known as FRAC3DVS. HGS is developed by Aquanty Inc. in Canada. 

In order to accomplish integrated analysis, HydroGeoSphere utilizes 
a rigorous, mass conservative modelling approach that fully couples the 
surface flow and transport equations with the three dimensional, variably 
saturated subsurface flow and transport equations. This approach is sig- 
nificantly more robust than previous conjunctive approaches that relied 
on the linkage of separate surface and subsurface modelling codes. HGS 
uses a globally implicit approach to simultaneously solve two-dimension- 
al diffusive wave equations and the three-dimensional form of Richards’ 
equation based on unstructured finite element grids. For each timestep, 
the model solves surface and subsurface flow, and solute and energy trans- 
port equations simultaneously, and provides a complete water and solute 
balance. HGS has the following features for surface and subsurface water 
modelling: 


e surface domain represented as two-dimensional 
overland flow; 


e subsurface domain consisting of three-dimensional 
unsaturated/saturated flow; 


e surface/subsurface domains interacting through 
physically based fluid exchange; 


e temporally and spatially varying evapotranspiration 
based on land use; 
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e impact of snowmelt on hydrologic regime; 
e delineation and tracking of the water table position; 


e handling of non-ponding or prescribed ponding 
recharge conditions and seepage faces; 


e representation of fractured geologic materials with 
arbitrary combinations of porous, discretely fractured, 
dual-porosity, and dual-permeability media for the 
subsurface; 


e accommodation of storage, solute mixing and variable 
flow distribution along wellbores; and 


e density-dependent flow and transport. 


The capabilities of the mass and heat transfer module of HGS include: 


e capability of modelling non-reactive and reactive 
chemical species transport in the associated surface 
and subsurface flow fields, including solute interactions 
between the surface and subsurface flow regimes; 


e calculation of temperatures in the surface and 
subsurface flow regimes as driven by air temperature 
and incoming solar radiation, accounting for land 
surface-atmospheric thermal interactions; 


e handling of fluid and mass/thermal energy exchanges 
between fractures and matrices, including matrix 
diffusion effects and solute/thermal energy advection in 
the matrix; and 


e straight or branching decay chains representing 
degradation reactions. 


HGS is written in FORTRAN and is being continuously developed. It 
runs on all versions of Windows and Linux systems. HGS does not cur- 
rently have a graphical user interface (GUI). All model parameters, grid 
structures, material properties, or numerical parameters are written in 
text files. A preprocessor (called GROK) prepares the input files for HGS. 
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The main input file is an instruction-driven text file that only requires a 
text editor. HGS has a post-processing program called HSPLOT to convert 
the output data to a format that can be read by third-party visualization 
packages such as TECPLOT or GMS (Aquanty, 2015). HGS has been used 
for simulating variably saturated groundwater flow and reactive solute 
transport, such as nitrate (Koh, Lee, & Lee, 2016). 


FEFLOW 

The Finite Element subsurface FLOW system (FEFLOW) is a two-dimen- 
sional/three-dimensional finite-element model for simulating ground- 
water flow, solute, and heat transfer in porous media and fractured media 
(DHI, 2016). 

FEFLOW uses finite-element-based analysis to solve the groundwater 
flow equation of both saturated and unsaturated conditions as well as 
reactive multi-species solute and heat transport, including fluid density 
effects and chemical kinetics for multi-component reaction systems. 
Contaminant transport processes include advection, hydrodynamic dis- 
persion, linear and nonlinear sorption isotherms, and first-order chem- 
ical, non-equilibrium, kinetic reactions between species. FEFLOW is 
available for Windows systems as well as for different Linux distributions. 
FEFLOW is a completely integrated package from simulation engine to 
user interface. The option to use and develop user-specific plug-ins via the 
programming interface (Interface Manager IFM) allows for the addition 
of external code or even external programs to FEFLOW. 

Since its birth in 1979, FEFLOW has been continuously extended and 
improved. It is consistently maintained and further developed by a team of 
experts at DHI-WASY. FEFLOW is used worldwide as a high-end ground- 
water modelling tool at universities, research institutes, government 
agencies, and consulting companies. FEFLOW can be efficiently used to 
describe the spatial and temporal distribution and reactions of ground- 
water contaminants (Regnery et al., 2017), to model geothermal processes, 
to estimate the duration and travel times of chemical species in aquifers, 
to plan and design remediation strategies and interception techniques, 
and to assist in designing alternatives and effective monitoring schemes. 
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Table 8. Summary of Selected Groundwater Quality Models 


Model & Spatial Type of Simulated Simulated GUI Availability References 
Source Representation Simulation Processes Parameters and Support 
MT3D-USGS 3D finite Steady-state Advection, Reactive and Yes Public Bedekar et al., 
(USGS) difference grid or transient dispersion, and non-reactive (public domain 2016; Zheng & 
solution chemical reactions contaminants domain Wang, 1999 
(sorption, first- and 
order decay, propri- 
interspecies ety) 
reactions, and 
parent-daughter 
chain reactions) 
of contaminants 
in groundwater 
systems coupled 
with MODFLOW. 
Transport within 
streams and 
lakes, including 
solute exchange 
with connected 
groundwater. 
SUTRA (USGS) 2D/3D finite Steady-state Saturated- Reactive and Yes Public Tsanis, 2006; C. 
element grid or transient unsaturated, non-reactive,  (par- domain Voss, 1999; C. |. 
solution fluid-density- single tially Voss & Provost, 
dependent ground- species publicly 2010 
water flow with avail- 
energy transport able, 
or chemically- partially 
reactive (sorption, propri- 
first-order decay) ety) 


single-species 
solute transport. 
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Table 8. (continued) 


3D finite element Steady-state 
grid or transient 
solution 


HydroGeo- 
Sphere 
(Aquanty) 


Entire terrestrial 
portion of the 
hydrologic cycle; 
uses a globally- 
implicit approach 
to simultaneously 
solve the 2D 
diffusive-wave 
equation and 

the 3D form of 


Richards’ equation. 


Reactiveand No Propriety 
non-reactive 


contaminants 


Aquanty, 2015; 
Koh et al., 2016 


MIKE SHE Square finite 
(DHI) difference grids: 
2D overland, 
1D channels, 
1D-unsaturated, 
and 3D 
saturated zones 


Steady-state 
or transient 
solution 


Entire hydrologic 
cycle, including 
groundwater and 
solute transport 
with simplified 
reaction processes 
(desorption, 
degradation, 

plant uptake, 
multi-species 
kinetics). Finite 
difference method 
for saturated 
groundwater 

flow, several 
representations of 
unsaturated flow, 
including the 1D 
Richards equation, 
MIKE 11 for flow in 
river and stream 
networks, and 

the 2D diffusive- 
wave approach for 
overland flow. 


Reactiveand Yes Propriety 
non-reactive 


contaminants 


DHI, 2017c, 
2017d; Refsgaard, 
Thorsen, Jensen, 
Kleeschulte, & 
Hansen, 1999 


PL 
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Table 8. (continued) 


Model & Spatial Type of Simulated Simulated GUI Availability References 
Source Representation Simulation Processes Parameters and Support 
FEFLOW (DHI) 2D/3D finite Steady-state Saturated and Reactive and Yes Propriety DHI, 2016; 


element grid 


or transient 
solution 


unsaturated 
density dependent 
flow, transport 

of mass (multiple 


solutes, desorption, 


degradation, 
kinetic reactions 
between species), 
and heat 


non-reactive 
contaminants 


Regnery et al., 
2017 


3.4 Land use/land cover models 


Land use/land cover (LULC) modelling helps to explain and/or predict 
LULC change processes (Pontius & Schneider, 2001) to understand the 
linkage between socioeconomic processes associated with land develop- 
ment and natural resource policies (Brown et al., 2000) and the causes 
and consequences of changes in the spatial and temporal patterns of land 
conversion (Irwin & Geoghegan, 2001). Furthermore, the models repre- 
sent a simplification of the complex behavior of the socioeconomic and 
physical environments, and ecological changes. Such models are used to 
explore future land-use changes under different scenarios and conditions 
(Veldkamp & Verburg, 2004). According to Lambin et al. (2000), LULC 
change modelling can address at least one of the following: 


1. Socioeconomic and environmental variables that mostly 
explain land-cover changes. 


2. ‘The locations that are affected by land-cover change. 


3. The rate at which land-cover changes progress. 


In terms of modelling, LULC deals with a complex structure of linkages 
and feedbacks to analyze the dynamics of LULC practices in the past, with 
the intention of determining trajectories of change and projecting possible 
future changes. 

The first generation of LULC models dates from the 1940s. The models 
considered an entire land system as a static entity. Land uses and activities 
were simulated at a cross-section in time, and their dynamics were con- 
sidered as trending toward self-equilibrium. These models were criticized 
because they had no spatial structure (Silva & Wu, 2012). Later on, the 
models started to take the spatial dimension into account with cellular 
automata (CA) models (Fredkin, 1990), and also to use GIS for data inte- 
gration and spatial analysis, such as the California Urban Futures Model 
(CUF) (Landis, 1995). 

Many LULC change models have been developed for different appli- 
cations, such as deforestation (Kaimowitz & Angelsen, 1998), agricultural 
intensification (Lambin et al., 2000), land-use based on economic theory 
(Bockstael & Irwin, 2000), and urban studies (Mitasova & Mitas, 1998). 
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Table 9. Land Use/Land Cover (LULC) Change Models 


Description 


Watershed studies 


Study area 


Model Reference 

Land-use change Ott & 

scenario kit (Luck) | Uhlenbrook, 
2004 


This model executes LULC change scenarios 
based on the characteristics of each grid cell 
and its relationships to neighboring cells. It 
has three modules for simulation of urbaniza- 
tion, agriculture, and forest LULC change. 
Luck includes the topography, soil, river 
network, and axes of infrastructural 
development. The spatially averaged, large- 
scale trend of LULC development must be 
provided as an external input, the so-called 
“scenario target.” 


Quantified the impact of LULC 
changes at the event and sea- 
sonal time scale. 


Dreisam basin in 
southwest Germany 


Land Transforma- Tang et al., 2005 


tion Model (LTM) 


A LULC forecasting model that employs a set 
of spatial interaction rules and machine learn- 
ing using artificial neural network to identify 
the nature of spatial interactions of driving 
forces based on the historical data to forecast 
future LULC change. 


Assessed the impacts of future 
LULC change on long-term 
runoff and non-point source 
pollution. 


Muskegon River 
watershed, in the 
eastern coast of Lake 
Michigan 


Conversion of Lin et al., 2007 
Land Use and its 


Effects (CLUE) 


This model is based on the spatial alloca- 
tion of demands for different LULC types to 
individual grid cells. It combines biophysical 
and human LULC drivers in space and time 
(Veldkamp and Fresco, 1996). The model is 
an empirical model which begins with the 
evaluation of relationships between land use 
and its driving factors and continues with the 
addition of dynamic simulation of inter- 
actions between the spatial and temporal 
dynamics of land use systems. 


Used the CLUE-s model (which 
is a modified version of CLUE) 
to simulate various future 

land use scenarios based on 
driving factors with spatial 

and non-spatial policies to 
assess the impact of future 
LULC change on hydrological 
processes. 


Wu-Tu watershed in 
northern Taiwan 


Įuəuodwop Ye 104 SOYDeOIddY BulljapoW O'S 


ZL 


Table 9. (continued) 


Urban growth Lin et al., 2008 
model Slope, Land 

use, Excluded 

land, Urban extent, 

Transportation, 

and Hillshading 


(SLEUTH) 


The model is a cellular automata class model, 
which is a probabilistic model that uses 
Monte Carlo routines to generate multiple 
simulations of urban growth. 


Applied SLEUTH and CLUE-s 
models to to analyze the effects 
of future urban sprawl on the 
LULC patterns and hydrological 
processes. 


Paochiao watershed 
in Taipei County, 
Taiwan 


Urban 
Development 
Simulation Model 


Cuo et al., 2011 


UrbanSim incorporates the interactions be- 
tween land use, transportation, the economy, 
and the environment. 


A land cover change mod- 
el (LCCM), UrbanSim, and 
biophysical site and landscape 


Puget Sound basin, 
Washington 


(UrbanSim) characteristics were used to 
investigate the potential im- 
pacts of projected future land 
cover and climate change on the 
hydrology. 
“What if?” McColl & Aggett, An easy-to-use GIS-based planning support ntegrated a land-use/cover Kittitas County, 
2007 system that can be used to explore the orecasting model with an Washington 
most important and difficult aspects of the event scale, rainfall-runoff 
land planning process: conducting a land model to improve LULC policy 
suitability analysis, projecting future land ormulation in the study area. 
use demand, and allocating the projected 
demand to suitable locations. 
Cellular Wijesekara etal., CA is a rigorous modelling approach for nvestigated the impact of Elbow River 


Automata (CA) 2012 


characterizing complex spatial systems 
through a bottom-up simulation of local 
interactions between neighboring cells. 


uture (20 years) land-use 
changes on the hydrological 
processes using a LULC cellular 
automata (CA) model and the 
distributed physically-based 
MIKE SHE/MIKE 11 hydrological 
model. 


watershed in southern 
Alberta 


NERC/ESRC Land Dunn & Mackay, 
Use 1995 
Programme 

(NELUP) 


The model applies a general system 
framework for organizing the large amounts 
of information that are relevant to decision- 
making in land use. 


Evaluated the potential 
impacts of land use change on 
evapotranspiration. 


Tyne River basin, UK 


However, little attention was paid to simulate future land-use/cover chan- 
ges for hydrological studies at a watershed scale. Most of them have con- 
sidered historical LULC maps as static input for land-related parameters 
in their simulations. In the recent past, very few studies have reported 
developing LULC modelling to forecast future LULC to analyze the effect 
of LULC changes on the catchment water balance (Wijesekara et. al., 2011; 
Cuo et al., 2011; Lin et al., 2007; Tang et al., 2005). Table 9 lists the widely 
used LULC models in the literature. 

In this report, we classify LULC models into two broad groups: empir- 
ical models and simulation integrated models. 


3.4.1 Empirical models 

Empirical models consider the statistical/mathematical relationship be- 
tween LULC classes, and external driving factors, such as physical char- 
acteristics (biodiversity, soil functions, and water resources), population, 
and technological development. However, they do not consider the inter- 
actions among LULC classes over time and the human behavior that 
can lead to the spatial process/outcome of land-use changes. The models 
can be implemented based on different statistical methods. For example, 
multivariate regression models (Braimoh & Vlek, 2004) and logit regres- 
sion models (Turner et al., 2001) have been used to evaluate the possible 
exogenous contributions of causal factors (Yu et al., 2011). They can also 
use some assumptions of possible future land developments to project fu- 
ture land-use patterns such as ‘what if? scenarios (Braimoh & Vlek, 2004). 

Empirical models can be spatial or spatially explicit models (Yu et al., 
2011). Spatial models have been primarily developed by economists to de- 
scribe spatial patterns of land use (Yu et al., 2011). These models are mostly 
based on the concept of distance to driving factors, and they ignore other 
important features of the landscape. 

According to Agarwal et al. (2002), there are two groups of spa- 
tially explicit models: spatially representative and spatially interactive. 
However, only spatially representative models belong to empirical models 
in our classification. These models deal with data in two or three spatial 
dimensions (for example, northing, easting, and elevation). These models 
cannot examine topological interactions between geographical features at 
each of the timesteps. However, the value of each feature can change or 
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stay constant independent of neighboring features. (Spatially interactive 
models, which incorporate spatial relationships and interactions between 
neighboring units over time, belong to the simulation integrated models 
class described in the next section.) 


3.4.2 Simulation integrated models 

Simulation integrated models are capable of conducting integrated simu- 
lations of LULC dynamics of landscape influenced by several factors. 
Integrated models represent “the relationships, interactions, and feed- 
backs between spatial and non-spatial components of a LULC system, 
such as human and economic activities, and physical and environmental 
characteristics of the landscape. On the other hand, simulation models 
are mathematical models that use computational resources to simulate 
the dynamics of land surface. However, in terms of modelling technique, 
Wilson (1974) explains that they are: 


a set of rules which enable a set of numbers to be operated 
upon, usually in the computer, although the rules and the 
consequences of applying them cannot be written down as 
a set of algebraic equations. . . . Sometimes, the simulation 
technique lends itself naturally to a problem. This happens, 
for example, when the underlying theory consists of a set 
of statements involving conditional probabilities. ... We 
resort to simulation techniques for situations which are 
too complicated to be handled by more straightforward 
algebraic techniques. 


The meaning of simulation methods of modelling was clarified by Batty 
(1976) as: 


Analytic methods of modelling use mathematical analysis 
to reach at explicit equations representing the behavior of 
the system while simulation methods are used to derive the 
behavior of the system when the system is too complex to be 
modeled using the more direct analytic approach. 
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In general, simulation integrated models derive the behavior of a complex 
system, and simulate the change of the land’s attributes using a set of rules 
based on the interactions, relationships, and linkages arising between 
components of land use either internally (the number of LULC classes 
and their interactions among neighborhood units) or externally, with the 
overall intent to capture the dynamics of land and reproduce its patterns 
over time. 


3.5 Climate models 


Climate models are important for improving our understanding and abil- 
ity to predict climate behavior as a result of natural variability and change 
and human activity, as well as understanding the impacts of climate change 
and variability on environmental processes (Farjad et al., 2019). Climate 
models are mathematical methods and computer programs which simu- 
late interactions between land, and/or atmosphere, and/or ocean, and/or 
ice by incorporating physical system processes. They are run using power- 
ful computers to capture the complex influence of internal processes and/ 
or external forcing on the climate system. They range from simple energy 
balance models to complex Earth Systems Models (ESMs). To keep it sim- 
ple, climate models could be divided into three categories: simple models, 
general circulation models, and intermediate complexity models. 


3.5.1 Simple models 

These are the earliest and most simplified version of climate models de- 
veloped based on the concept of energy balance. The simple linear relaxa- 
tion model (usually referred to as an Energy Balance Model) is the simplest 
climate model in this category and is time dependent. These models can 
range from zero- to two-dimensional models. Zero-dimensional models 
assume a balance between incoming solar radiation and outgoing long- 
wave radiation, resulting in a uniform temperature over the land surface. 
One-dimensional models assume different latitudinal zones cover land 
surfaces with different incoming solar energies. However, two-dimension- 
al models assume variations in two directions. 
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3.5.2 General circulation models 
General circulation models, also known as Global Climate Models 
(GCMs), are numerical models used to simulate the complex interaction 
of the processes and feedbacks in the climate (M. Wang et al., 2012). They 
are developed based on physical laws of climate dynamics, and they are 
solved by mathematical equations or sometimes empirical relations. For 
example, atmospheric GCMs numerically solve the equations of phys- 
ics, such as dynamics, radiative transfer, or thermodynamics, as well as 
chemistry applied to the atmosphere and its constituent components. 
Whereas in more primitive GCMs only the thermodynamic role of the 
ocean was considered, GCMs today typically include the dynamics of the 
ocean and its interactions with the atmosphere and are therefore known 
as Atmosphere-Ocean GCM (AOGCM) or coupled atmosphere-ocean 
models. The term Global Climate Models (GCMs) is typically used to refer 
to climate models that reflect both atmospheric and oceanic processes 
and feedbacks. The current generation of GCMs includes the hydrologic- 
al cycle (which couples terrestrial, atmospheric, and ocean reservoirs of 
water and the flows between these reservoirs), terrestrial biosphere, con- 
tinental ice sheets, and the ocean's carbon cycle and its interactions with 
the atmosphere and the ocean. As a result, a variety of climate compon- 
ents are included such as atmospheric and ocean circulation, atmospheric 
temperature profiles, snow and ice distribution, and wind patterns. 
General circulation models are numerical models used to simulate the 
complex interaction of the processes and feedbacks in the atmosphere (M. 
Wang et al., 2012). They are developed based on physical laws of climate dy- 
namics of the atmosphere, and they are solved by mathematical equations 
or sometimes empirical relations. Primitive GCMs prescribed the physical 
characteristics of the earth’s surface (e.g., sea-surface temperature, land 
temperature, and soil wetness) as boundary conditions. Subsequently, 
coupled atmosphere-ocean models have been developed to reflect inter- 
actions between the atmosphere and the ocean. The term Global Climate 
Models (GCMs) is typically used to refer to climate models that reflect both 
atmospheric and oceanic processes and feedbacks. GCMs today typically 
include the dynamics of the ocean and its interactions with the atmos- 
phere and are therefore known as Atmosphere-Ocean GCM (AOGCM) 
or coupled atmosphere-ocean models. For example, atmospheric GCMs 
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numerically solve the equations of physics, such as dynamics, radiative 
transfer, or thermodynamics, as well as chemistry applied to the atmos- 
phere and its constituent components. The current generation of GCMs 
includes the hydrological cycle (which couples terrestrial, atmospheric, 
and ocean reservoirs of water and the flows between these reservoirs), ter- 
restrial biosphere, continental ice sheets, and the ocean’s carbon cycle and 
its interactions with the atmosphere and the ocean. As a result, a variety 
of climate components are included such as atmospheric and ocean circu- 
lation, atmospheric temperature profiles, snow and ice distribution, and 
wind patterns. 

Further, a subset of climate modelling involves Integrated Assessment 
Models (IAMs) by incorporating socioeconomic aspects to understand 
how societal factors influence the climate (for example, how population, 
economic growth, and use of fossil energy influence the climate on earth). 
Therefore, IAMs produce scenarios of future greenhouse gas emissions, 
and these scenarios are then run through various Earth System Models 
(ESMs) to generate future climate projections, used for assessing the im- 
pact of climate change to develop adaptation and mitigation strategies and 
policies for sustainable growth and environmental stewardship for future 
generations. 

GCMs provide coarsely scaled outputs in spatial and temporal reso- 
lutions, as they are computationally intensive. Therefore, these scales are 
nearly inadequate for direct use in environmental models (S. Tripathi et 
al., 2006; Sunyer et al., 2012), and so the outputs of a low-resolution climate 
model need to be downscaled to a finer suitable scale (Teng et al., 2012), 
which can be dynamical or statistical (Sunyer et al., 2012). Downscaling 
methods are an important tool in assessing the impact of future climate 
change on environmental processes at both regional and local scales. In 
fact, downscaling methods bridge the gap between the resolution of climate 
models and environmental models (Fowler et al., 2007). Various studies 
have used downscaling methods to produce the required meteorological 
variables for environmental modelling (X.-C. Zhang, 2005; Chen et al., 
2010; Willems & Vrac, 2011; Li et al., 2012; Farjad et al., 2015). Dynamical 
downscaling refers to the use of process-based regional climate models 
(RCMs) to provide climate data within boundary conditions prescribed by 
a GCM, through regional-scale atmospheric simulations, at a finer spatial 
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and temporal resolution than GCMs (Phatak et al., 2011). On the other 
hand, statistical downscaling relies on the statistical relationships between 
large-scale climate model variables (GCMs or RCMs) and local-scale cli- 
mate variables (Sunyer et al., 2012). 


3.5.3 Intermediate complexity models 

Intermediate complexity models strike a balance between simple energy 
balance models and GCMs, in terms of degree of complexity. These models 
are also known as Earth Models of intermediate Complexity (EMICs). In 
terms of dynamics and resolution, these models are simpler than GCMs, 
but they are more comprehensive than simple models, in terms of the 
number of components and processes (e.g., LOVECLIM). Some of these 
models are built for specific purposes with a specified range of atmos- 
pheric components, but most of them include sea ice, dynamic vegetation, 
land surface processes, and ice sheet models. LOVECLIM1.2 incorporates 
atmosphere, land surface, ocean and sea ice, ice sheets/icebergs, and the 
carbon cycle. 


3.6 Ecological models 


Ecological models are mathematical models for biological and biophysic- 
al processes (which can be analytic or simulation-based) to understand 
ecological processes and capture changes in ecosystems. Most ecologic- 
al models are generally integral parts of core system analyses or water- 
shed models. For example, the AQUATOX model is an integral part of 
the BASINS system with links to the watershed models SWAT and HSPF, 
whereas ECO Lab is an integral part of MIKE SHE/MIKE Hydro. The 
widely used ecological models are described as follows: 


3.6.1 AQUATOX model 

The AQUATOX model simulates the fate of sediments, organic chemicals, 
and nutrients, and their influence on the ecosystem (such as fish, inver- 
tebrates, or aquatic plants) and thus is capable of integrated modelling of 
ecology/biology and water quality. For example, Bingli et al., (2008) used 
the AQUATOX model to simulate the environmental fate and aquatic 
ecological impacts of a nitrobenzene concentration in the Songhua River, 
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China. They conducted a sensitivity analysis to determine the key process- 
es that influence the nitrobenzene concentration levels and found signifi- 
cant changes in biomass for diatoms and mussels. Other typical applica- 
tions of the AQUATOX model are: 


e calculating recovery time of contaminated fish tissues to 
safe levels when pollutant loads are reduced, 


e evaluating impacts of pesticides and other toxic 
substances, 


e developing numeric nutrient targets based on desired 
biological endpoints, 


e evaluating effects of multi-stressors on biological 
systems, and 


e measuring the impact of climate change on the 
ecosystem. 


3.6.2 Mike ECO Lab 

MIKE ECO Lab is a generic ecological modelling tool for simulating pro- 
cesses related to water quality and ecological systems. One of the advan- 
tages of the MIKE ECO Lab compared to other ecological models is that 
it can be linked to the range of one-dimensional, two-dimensional, and 
three-dimensional MIKE modelling systems to address a variety of com- 
plex ecological issues. In addition, MIKE ECO Lab not only contains a 
generic equation solver, but it can also be applied as a generic post-proces- 
sor of hydrodynamic results, for instance, calculating flood risk indices or 
a scour risk formula. Santos et al. (2015) used the MIKE ECO Lab linked 
with Mike Hydro Basin model for water quality assessment, such as the 
dissolved concentration of phosphorus, in a river basin with recurrent 
wildfires in northern Portugal. They found a positive correlation between 
the occurrence of forest fires and the concentration of phosphorus in the 
water. Other typical applications of the MIKE ECO Lab model include: 


e study of simple and complex ecological systems, 


e water quality modelling related to surface/subsurface, 
rivers, wetlands, lakes, reservoirs, estuaries, coastal 
waters, and the sea, 
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e modelling of ecosystem response spatially, and 


e impact and remediation assessment. 


3.6.3 BASS 

Bioaccumulation and Aquatic System Simulator (BASS) is used to simulate 
the bioaccumulation of chemical pollutants and dynamics, and the popu- 
lation of fish assemblages under chemical and non-chemical stressors. 
BASS is a process-based model that can simulate ecological, physiological, 
and toxicokinetic processes, and can address the limitations of simple bio- 
accumulation factor (BAF) approaches for predicting concentrations of 
extremely hydrophobic chemicals and metals. Knightes et al. (2009) used 
BASS to model fish mercury-intake as a function of gill exchange and 
dietary ingestion. The model partitions mercury internally to water, lipid, 
and non-lipid organic materials. Knightes et al. estimated a physiologic- 
ally based carrying capacity for zooplankton and phytoplankton based on 
projected oxygen consumption and prevailing dissolved oxygen content. 
Typical applications of the BASS model are: 


e simulating fish methylmercury bioaccumulation, 


e estimating lag times of mercury residues in fish in 
response to mercury load reductions, 


e simulating time dynamic bioaccumulation when simple 
steady-state methods (e.g., BSAFs or BAFs) are not 
considered sufficient, and 


e evaluating responses of fish community composition, 
production, biomass, and PCB/mercury bioaccumulation 
potential to changes in climate, LULC, and fisheries 
management scenarios. 


3.6.4 SERAFM 

SERAFM is a process-based (steady-state) mercury cycling model used to 
estimate mercury concentrations in the water column, fish tissue, and sedi- 
ment for the species Hg0, HgII, and MeHg. Sub-modules of SERAFM con- 
sist of mercury loading (watershed and atmospheric deposition), abiotic 
and biotic solids balance (soil erosion, settling, burial, and resuspension), 
equilibrium partitioning, water body mercury processes, and wildlife risk 
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calculations. S. Brown et al. (2007) used the SERAFM model to estimate 
unfiltered and dissolved total mercury concentration, as well as MeHg, in 
the water column in Steamboat Creek in Nevada. They highlighted the 
capability of the SERAFM model of estimating mercury concentration in 
arid western environments. Typical applications of the SERAFM model 
are: 


e mercury deposition impacts on a water body, 


e estimating steady-state mercury cycling in a river, lake, 
and watershed, 


e wildlife risk prediction exposed to mercury- 
contaminated sediments, and 


e evaluating sensitivity of model parameters and processes 
to pollution. 


3.6.5 Physical Habitat Simulation System (PHABSIM) 

The Physical Habitat Simulation System (PHABSIM) model is used to 
simulate relationships between river flow and physical habitat for a var- 
iety of life stages of a species of fish or a recreational activity. PHABSIM 
integrates a river model with a biological model of habitat (based on habi- 
tat suitability criteria) to estimate changes in a habitat index (weighted 
usable area) as a function of river discharge. Some of applications of the 
PHABSIM model are: 


e predicting the micro-habitat conditions in rivers and the 
relative suitability of those conditions to aquatic life, 


e understanding the impact of mining-derived sediment 
on aquatic physical habitat, 


e determining the utility of a reconnaissance-level 
physical habitat suitability, and 


e examining the trade-off between the value of water used 
instream with the water used out-of-stream. 
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3.7 Air quality models 


Understanding the impact of emitted air pollutants from natural or man- 
made sources is a challenging topic in CEA studies. 

There are many approaches to quantify air pollution, and they can be 
categorized as dispersion, photochemical, and receptor models, as well as 
geospatial-based models. Different factors are considered in the literature 
for selecting a proper approach, such as availability of data, temporal and 
spatial scale, acceptable level of uncertainty, type of the pollutant, and type 
of the activity (biomass burning, traffic-based pollution). Here, the most 
commonly used approaches to investigate air pollution are explained. 


3.7.1 Dispersion modelling 

Dispersion models are mathematical models that predict the concentration 
of pollutants at specified ground level receptor locations. These models use 
emissions and meteorological variables as inputs. The most common dis- 
persion models are AERMOD and CALPUFF. The AERMOD model in- 
corporates air dispersion based on a planetary boundary layer turbulence 
structure and scaling concepts while, the CALPUFF model simulates the 
effects of space- and time-varying meteorological conditions on pollution 
transport, transformation, and removal. Other common dispersion mod- 
els are CALINE3, BLP, CAL3QHC/CAL3QHCR, OCD, and CTDMPLUS. 


3.7.2 Photochemical modelling 

Photochemical models simulate atmosphere pollutant concentrations 
using a set of mathematical equations characterizing the chemical and 
physical processes over large spatial scales. The most common photochem- 
ical models are the Community Multiscale Air Quality model (CMAQ), 
the Comprehensive Air quality Model with extensions (CAMx), and the 
Regional Modelling System for Aerosols and Deposition (REMSAD). 


3.7.3 Receptor modelling 

Receptor models are mathematical or statistical procedures that use the 
measured physical and chemical characteristics of gases/particles (not 
pollutant emissions and meteorological data, unlike photochemical and 
dispersion models) to identify and quantify the source of pollutants at 
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receptor locations. The most common receptor models are the Chemical 
Mass Balance (CMB) model, Unmix model, and Positive Matrix 
Factorization (PMF) model. 


3.7.4 Geospatial-based modelling 
These models typically generate the spatial distribution of pollutants and 
can be categorized as: 


3.7.4.1 Proximity-based models: 

These models are GIS based and are used for assessment of exposure to air 
pollution when the density of spatial monitoring stations is sparse. These 
models use simple approximate measurements such as ‘buffering’ to as- 
sess the exposure, which results in severe limitations on their usefulness. 


3.7.4.2 Interpolation-based models: 

Spatial interpolation models have been developed to estimate values at 
unknown locations based on known values (i.e., measurements) by using 
deterministic and stochastic geostatistical techniques. The interpolation 
of values between monitor locations is an approach that can be completed 
entirely within GIS. These models use different geostatistical interpola- 
tion techniques. One conventional technique is inverse distance-weight- 
ed interpolation (IDW). For any given location in the study area (i.e., a 
residential address or postal code), a weighted average concentration is 
developed, with the highest weight given to the nearest monitors, there- 
by producing a continuous pollution surface. These methods assume a 
stronger correlation among points that are close together versus those that 
are farther apart. 

The most popular interpolation-based model used in air pollution 
studies is Kriging, an optimal interpolation technique that makes the best 
linear unbiased estimate of the variable’s value. Interpolation models use 
pollution measurements, which offer primary advantages over the prox- 
imity models. However, the main disadvantages of interpolation models 
are that they depend on the availability of monitoring data and require a 
dense network of sampling data. 
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3.7.4.3 Land Use Regression (LUR) models: 

Land Use Regression (LUR) models are multivariate regression models 
that are used for estimating individual exposure to air pollution at a fine 
spatial scale. They can estimate the pollution concentration at any given 
location by using surrounding attributes such as LULC classes, traffic, and 
topography within the area as independent variables (x). Therefore, the 
measured levels of pollutants in LUR models are considered as dependent 
variables (y). 
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4.0 INTEGRATED ENVIRONMENTAL 
MODELLING 


Integrated Environmental Modelling (IEM) allows the environment to 
be considered in a holistic way and provides a science-based system to 
explain (the past) and predict (the future) behavior of environmental sys- 
tems in response to human and natural sources of stressors. IEM often 
requires integrating (spatial) data and computational models from a var- 
iety of disciplines (e.g., related to physical, biotic, social, and economic 
environments) and at different scales, to understand and to solve complex 
societal problems that arise from the interaction of humans and the en- 
vironment, and to contribute in this way to establishing the foundation of 
sustainable development, to inform policy, and to support decision-mak- 
ing (Rothman, 1997; Parker, 2002). 

Model integration is achieved by linking together stand-alone mod- 
els or model components using various coupling approaches. Coupling 
approaches are defined from various perspectives, such as the degree and 
direction of linkage. Standard classifications and names have not yet been 
established. There are a variety of coupling approaches used in this review, 
which are classified based on the calculation order of model components: 


e Fully coupling: Equations governing all model 
components are solved simultaneously within a single 
monolithic code. 


e Dynamic coupling: Two or more individual models are 
tightly coupled via the exchange of data dynamically 
during simulation at each timestep/predefined 
frequency. 
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e Sequential coupling: Two or more individual models are 
sequentially run and loosely coupled via the modelling 
output files by which the output of one model forms the 
input of the other. 


e Interactive coupling: Two or more individual models are 
run in an alternative order following predefined time 
periods and loosely coupled via the modelling output 
files by which the output of one model forms the input of 
the other. 


e Iterative coupling: Two or more individual models are 
loosely coupled by which one model is iteratively called 
as a slave to provide the output of simulation results to 
feed another master model. 


e Hybrid coupling: Component models are integrated with 
two or more coupling approaches. 


The key advantage of fully coupled models is there capacity to solve compo- 
nent models with concurrent feedbacks from each other without a delay in 
timesteps. There are fewer examples of fully coupled models at regional or 
watershed scales. The applications mentioned in literature are mainly for 
integrated subsurface and surface flow and solute transport models, such 
as MIKE SHE (Farjad et al., 2017a; Farjad et al., 2017b), HydroGeoSphere 
(Therrien et al., 2010), OpenGeoSys (Kolditz et al., 2012), and Parflow 
(Kollet & Maxwell, 2006). A dynamically coupled model transfers infor- 
mation at each simulation timestep and acts as a single model or frame- 
work component that communicates through feedback mechanisms, such 
as the MIKE SHE coupled with other models (e.g., MIKE 11) or a set of 
add-ons, including with DAISY and ECO Lab. Conversely, loosely coupled 
models lack feedback mechanisms during runs of component models and 
exchange information through file transfer mechanisms at or after each 
run of components. 

Many stand-alone deterministic hydrological and water quality mod- 
els are available. However, matching model capabilities with the complex- 
ities of natural and engineered systems is a challenge. In recent years, the 
difficulties in transferring information between models prompted the 
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development of integrated modelling tools. General computer frameworks 
are available as integration tools for coupling component models, such 
as the Open Modelling Interface (OpenMI) developed by a consortium 
of European universities and private companies (Gregersen, Gijsbers, & 
Westen, 2007), the Object Modelling System (OMS) developed by the United 
States Department of Agriculture (Ahuja, Ascough II, & David, 2005), 
and FRAMES (Framework for Risk Analysis Multi-Media Environmental 
Systems) developed by the USEPA (Babendreier & Castleton, 2005). Figure 
7 illustrates an ecological modelling system for assessing the impacts of 
multiple stressors on stream and riverine ecosystem services within river 
basins, utilizing FRAMES to combine component models. Although sur- 
face and subsurface flow systems are naturally connected, they are often 
divided into separate compartments due to computational burdens as well 
as different temporal and spatial scales of the processes involved. This ap- 
plies to separation of flow between the unsaturated and saturated zone, 
as well as to surface (i.e., overland and channel) and groundwater flow. 
Although it sounds ideal to solve the different partial differential equa- 
tions of all component models simultaneously at each timestep in a fully 
coupled system, it may be computationally inefficient, as feedback from 
some components to others is relatively slow. It is reasonable to simulate 
interactions between the unsaturated zone and the groundwater for a 
watershed with a shallow groundwater table by using dynamic coupling 
at certain timestep intervals. For a watershed with a deep groundwater 
table, where the roots of plants and agricultural practices cannot reach, 
there is no physical reason to go beyond sequential, interactive, or iterative 
coupling. 
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Figure 7. Integrated Ecological Modelling System for the Coal River 
Basin (from Johnston et al., 2017). 


MLM 
SWAT HSI 
oa IW 
WASP \ pass 
ESP 


Note: SWAT = Soil Water Assessment Tool, MLM= Mercury Loading Model, WASP = Water Quality 
Analysis Simulation Program, HSI = Habitat Suitability Index, PiSCES = Piscine Stream Community 
Estimation System, BASS = Bioaccumulation and Aquatic System Simulator, ESP = Ecosystem 
Services Processor. 


4.1 Integrated surface water-groundwater 
quantity modelling 


Understanding complex interactions between surface and groundwater 
hydrology is challenging due to the nonlinear hydrodynamic nature 
of surface and subsurface components, particularly in heterogeneous 
conditions. This can be even more complex when it comes to the math- 
ematical representation of these interactions in a modelling system. 
There are a growing number of watershed models capable of simulating 
integrated surface and subsurface interactions, such as ParFlow, MIKE 
SHE, CATHY, HydroGeo-Sphere (HGS), PAWS, OpenGeoSys (OGS), 
PIHM, Cast3M, ATS, GEOtop, and tRIBS+VEGGIE. However, there are 
few intercomparison studies available in the literature involving regional 
scales (10° to 10° km’). For example, Maxwell et. al. (2014) and Kollet et. 
al. (2017) conducted an integrated watershed intercomparison study on 
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a series of benchmarking problems. Maxwell et al. (2014) compared sev- 
eral integrated watershed models: CATHY, HGS, OGS, PAWS, ParFlow, 
PIHM, and tRIBS+VEGGIE. The models simultaneously solved adapted 
forms of the Richards and shallow water equations based on three-di- 
mensional or mixed (one-dimensional vadose zone and two-dimension- 
al groundwater) formulations for subsurface flow and one-dimensional 
(rill flow) or two-dimensional (sheet flow) conceptualizations for surface 
routing. Kollet et al. (2017) used the same approach but a slightly differ- 
ent experiment using the MIKE SHE, ATS, CATHY, Cast3M, GEOtop, 
HGS, and ParFlow models. Overall, both studies found good agreement 
between models, especially for simple test cases, whereas some differences 
were identified that were mostly associated with mathematical and num- 
erical representation or in the parameterization of physical processes. This 
intercomparison might not be valid for regional scales due to the hetero- 
geneity of landscapes, topography, climatic, geomorphology, stream pat- 
terns, density, and geological units. 

The Ontario Ministry of Natural Resources (2011) conducted an in- 
tegrated watershed intercomparison study for regional scales of the most 
popular integrated models: GSFLOW, MIKE SHE, HydroGeoSphere, 
ParFlow, and MODHMS. The limitations of each model are summarized 
in the following section. 


GSFLOW 


e Empirical water budget formulation: While Precipitation- 
Runoff Modelling System (PRMS) includes a variety 
of methods for simulating surface water hydrologic 
processes, not all methods are enabled for integrating 
MODFLOW in GSFLOW. The GSFLOW implementation 
of PRMS represents the water interchange between the 
surface soil zone using three reservoirs: preferential 
flow, gravity flow, and capillary. The soil zone exchanges 
flow with the MODFLOW unsaturated zone, and 
the rate of interchange between these reservoirs is 
modelled empirically. However, identification of optimal 
parameters was found to be difficult when completing 
the case studies. 


4.0 Integrated Environmental Modelling 95 


e Restricted surface water time-stepping and hydraulic 
routing: The GSFLOW implementation of MODFLOW 
and PRMS does not allow for the timesteps in a surface 
water model to be less than one day. This limitation may 
influence the simulation of hydrologic processes such as 
runoff or infiltration and snowmelt, all of which occur 
during shorter periods (i.e., sub-daily) within a day. 
Also, the model cannot represent overland flow routing 
and complex hydraulic structures, which are important 
to properly represent surface water flow events that 
occur during short time periods. Although GSFLOW 
may be calibrated to account for longer-term hydrologic 
trends, it should not be considered suitable for many 
short-term events. 


MIKE SHE 


e Uniform Grid Resolution: The overall capabilities of 
MIKE SHE would be more advanced if a variable 
resolution grid system were present. This would allow 
grid refinement near features of importance such as wells 
and surface water bodies as well as regions of highly 
variable topography. From a computational perspective, 
this would also be beneficial as it would allow for more 
efficient application of computing resources (e.g., fine 
model resolution within areas of interest and coarse 
resolution in surrounding regions). 


e Source code: The source code is proprietary and not 
available for examination or modification. 


e Purchase price: The purchase price of the code is 
considered to be high as compared to other alternatives. 
However, the experience gained when completing the 
case studies demonstrated that the purchase price of the 
code can be offset on a single project by the time savings 
realized by having the user interface available, as well as 
by the overall flexibility offered by MIKE SHE. 
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HydroGeoSphere 


MODHMS 


ParFlow 


Computational effort: HydroGeoSphere’s simulation 
times may be on the order of weeks for a single scenario, 
which is not practical for many applications. 


Surface water hydrologic processes and features: 
HydroGeoSphere does not fully account for hydrologic 
processes such as snowmelt and hydraulic structures. 


Lack of a graphical user interface: While processing 
tools are available for components of the model (e.g., 
finite element mesh), there is not a single and complete 
graphical user interface available for HydroGeoSphere, 
which limits its ability to be cost-effective for most 
applications. 


Application in cold regions: The model does not include 
winter processes. Snowmelt is arguably one of the most 
important hydrological processes in cold regions. 


Source code: The source code is proprietary and not 
available to the public for examination or modification. 


Flexibility: The model is not flexible in terms of representing 
various hydrological processes at different levels of 
complexity (e.g., representing groundwater flow using a 
linear reservoir approach when subsurface data is sparse). 


ParFlow is primarily a research code that requires third- 
party software to visualize most of its output. 


ParFlow simulations cannot incorporate hydraulic 
structures (e.g., dams, weirs, etc.). 


Besides the above-mentioned comparison of these models, the MIKE 
SHE/MIKE 11 model offers two specific benefits which are useful for CEA: 
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G) MIKE SHE/MIKE 11 is capable of tracing water and pollutants from the 
ground surface, through the soil and groundwater, and back into the sur- 
face water, (ii) the MIKE package includes some other tools which can be 
linked to the integrated modelling system for a specific CEA application. 
For instance, MIKE ECO Lab and FEFLOW can be linked to the water- 
shed model for ecological modelling and advanced localized groundwater 
quality simulation, respectively. 

It should also be noted that while the MIKE SHE/MIKE 11 model is 
a useful tool for CEA in small to large scale (domain) regions, it might 
not be capable of supporting large regions when a high spatial resolution 
(grid size) configuration (e.g., < 200 m) is required. This is because, al- 
though MIKE SHE supports parallelized computations, the parallelized 
approach relies on a shared memory approach (OpenMP), which has lim- 
ited opportunities for scaling. MIKE SHE’s code can take advantage of 
multi-core processors, however, at a certain point (approximately eight 
multi-core processors) the cost of communication between parallel pro- 
cesses exceeds the benefits of additional processor cores. Furthermore, 
single model runs cannot be distributed across multiple computers. In 
this regard, high-resolution integrated models, called hyper-resolution 
models, have been developed. These models could enable more realistic 
process-level simulations that are critical for many important CEA ap- 
plications at high spatial and temporal resolutions. Hyper-resolution 
modelling requires large parallel-clustered computing resources and solu- 
tion algorithms that efficiently use these resources. While such systems 
are increasingly becoming available to scientists and modellers in many 
earth science disciplines (e.g., climate modelling) for continental/global 
scales, environmental modelling communities have been slow to utilize 
these computational resources for regional scales. Recently, a few studies 
(Kollet et al., 2010; Maxwell, 2013; Maxwell et al., 2015) have attempted to 
apply parallel and high-performance computing techniques to simulate 
surface water and groundwater interactions using high resolution mod- 
els. For example, Kollet et al. (2010) used an integrated model (ParFlow) 
to simulate the interactions between land surface processes and variably 
saturated flow in a heterogeneous subsurface for a maximum number of 
approximately 8x10° grid cells. The parallel performance of the model 
was investigated based on a scaling assessment on the JUGENE massively 
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parallel supercomputer. JUGENE is an IBM Blue-Gene supercomput- 
er with a total of 294,912 processors and 144TB of memory capable of 
0.825 PetaFLOPS (floating point operations per second) and is currently 
ranked the fourth fastest supercomputer in the world. They indicated that 
regional scale hydrologic simulations on the order of 103 km? are feasible 
at hydrologic resolution of ~10°-10' m laterally and 10°-10~ m vertically, 
with reasonable computation times, which had been previously assumed 
to be an intractable computational problem. The advantage of develop- 
ing high-resolution integrated cumulative predictive models is not only 
motivated by the potential of coupling different environmental processes 
for cumulative effects assessment, but also because these types of models 
may be useful in serving as virtual laboratories or realities. 


4.2 Integrated watershed and receiving water 
quality modelling 


Although watershed models may have both hydrologic and water qual- 
ity modules, their capacities to simulate hydrodynamic and water quality 
processes in receiving water bodies (e.g., rivers) are generally limited and 
often of one-dimensional and quasi-dynamic state. There is a need to take 
advantage of the capacities of existing receiving water models in simulat- 
ing complex hydrodynamics and water quality processes by combining 
watershed and receiving water models. In such an integrated model, the 
watershed sub-model simulates and provides discharges (i.e., surface run- 
off) and loadings (i.e., non-point) to the receiving water quality model, 
coupled with either a dynamic or a sequential linkage. Some applications 
are summarized in Table 10 as examples. 

For dynamic coupling, OpenMI is commonly reported as a tool for 
interfacing component models from different disciplines or domains. In 
the Pinios River catchment in Greece, approximately 10,500 km’, two al- 
ternative integrated models were developed (Makropoulos et al., 2010) for 
water quality evaluation using OpenMI. The first consisted of the rainfall 
runoff NAM module of MIKE 11, the hydrodynamic model RISH-1D, and 
the water quality model RISQ-1D, while the second used NAM, the MIKE 
11 hydrodynamic module, and the water quality model OTIS. The same 
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Table 10. Examples of Integrated Watershed and Receiving Water Quality Modelling Studies 


Coupling Author Region/watershed scale Watershed Receiving water Coupling tool & data 
application model quality model exchanges 
Dynamic Makropoulos et Pinios River catchment NAM MIKE 11 & OTIS; or, OpenMI, links for node 
al., 2010 (10,500 km?), Greece, RISH-1D & RISQ-1D connection, flow, water level, 
integrated hydrologic, BOD concentrations. 
hydraulic, and water quality 
modelling 
Shrestha et al., River Zenne, Belgium, SWAT SWMM OpenMI, SWAT output as 
2013 integrated sediment transport upstream boundary condition 
modelling for SWMM model. 
Mentzafou & Evros river basin, 2,778km?, MIKE SHE MIKE 11 Add-on in single code, 
Dimitriou, 2011 Greece coupling flow, recharge, 
nitrate concentrations. 
Malek- Upper East Fork Poplar Creek MIKE SHE MIKE 11+ ECOLAB Add-on in single code, 
Mohammadiet watershed, Tennessee coupling flow, recharge, TSS, 
al., 2012 and mercury concentrations. 
Sequential Michael Baker Illinois River watershed, HSPF EFDC (3D) User defined linkages 
Jr. Inc. et al., Oklahoma via output and input text 
2015 files. The HSPF model 
hourly results are used to 
provide streamflow, water 
temperature, suspended 
solids (TSS), organic carbon, 
nutrients (N, P), algae 
biomass, and dissolved 
oxygen as input data for the 
EFDC lake model. 
Sutula et al., Santa Margarita River HSPF EFDC + WASP (3D) User-defined linkages via 
2016 Watershed, California, output and input text files, 


nutrient management 


coupling hourly flow, nutrient 
loads. 
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Table 10. (continued) 


Huang et al., Ribble catchment, Northwest HSPF (rual RMN 1D (river) + User defined linkages via 
2017 England, 12,920 km? area) + EFDC (2D estuary) output and input text files, 
Infoworks coupling flow and E. coli. 
(urban area) 
+ DMHSF 
(for fecal 
indicator) 
Privette et al., Reedy River watershed, South LSPC WASP (3D) User defined linkages via 
2015 Carolina output and input text files, 
coupling hourly flow, total 
phosphorus, and total 
nitrogen. 
Shabani et al., Devils Lake watershed, North SWAT CE-QUAL-W2 User-defined linkages via 
2017 Dakota output and input text files, 
coupling daily flow and sulfate 
loads. 
Yue & Cobb Creek Watershed, SWAT EFDC User-defined linkages 
Derichsweiler, Oklahoma via output and input text 
2005 files, coupling daily flow, 
Chlorophyll-a, CBOD, Nitrate, 
Organic N, Mineral P, and 
Organic P loads. 
J.M. Johnston Albemarle-Pamlico SWAT +WMM WASP USEPA’s FRAMES was used to 
et al., 2011 Watershed, North Carolina (watershed define linkages, coupling daily 
and Virginia mercury flow, nutrients, and mercury. 
model) 
Mankin et al., Melvern Lake watershed, AGNPS EUTROMOD User-defined linkages via 


1999 


Kansas 


output and input text files, 
coupling annual flow and 
nutrient loads. 


pollution loads for both diffusive and point sources were assumed in both 
integrated models, and the same BOD decay coefficient and dispersion 
coefficient were used in both RISQ-1D and OTIS. The comparative analy- 
sis of the two configurations illustrated the significant differences for two 
river nodes between model components and (consequently) model results, 
even when OpenMI was used as the integrating medium and the model 
schemes were set up for the same study area by collaborating modelling 
teams. It is suggested that this variation is therefore a measure of the un- 
certainty related to the input data discrepancies and different modelling 
techniques. The visualisation of this significant uncertainty may be very 
important for decision-making, including but not restricted to the iden- 
tification of the required level of water treatment for local communities. 

An OpenMI-based integrated model was developed for the purpose 
of simulating the sediment dynamics for the River Zenne in Belgium 
using SWAT to model water and sediment fluxes from rural areas and 
SWMM to simulate the hydraulics of the river, canal, and sewer sys- 
tems in the downstream urban catchments (Shrestha et a., 2013). The 
SWAT model essentially formed the upstream boundary condition for 
the SWMM model. 

MIKE SHE is fully and dynamically integrated with a channel flow, 
transport code MIKE 11, water quality, and the ecological module ECO 
Lab. The exchange of surface and subsurface water and the loadings be- 
tween the two components take place during the whole simulation run. 
As the MIKE SHE/MIKE 11 system is a ready-to-use commercial package 
without a need for integration programming efforts, it can be directly ap- 
plied to simulate groundwater, surface water, sub-subsurface interactions, 
receiving water hydrodynamics, and advection/dispersion processes, 
while the water quality kinetics are simulated using ECO Lab in a single 
model. For example, the system has been applied to analyze the mercury 
cycle in the environment and provide forecasting capabilities for the fate 
and transport of contamination within the Upper East Fork Poplar Creek 
watershed in Tennessee (Malek-Mohammadi et al., 2012) and the trans- 
port and fate of nitrate in the Evros River basin in Greece in a large area 
about 2,778km? (Mentzafou & Dimitriou, 2011). 

In scientific and gray literature, many integrated watersheds and 
receiving water quality models developed using a sequential coupling 
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approach have been reported. For integrated water quality modelling, 
HSPF and SWAT are the two most selected watershed models, while 
MIKE 11, EFDC, and WASP are frequently selected hydrodynamic and 
water quality models (Huang et al., 2017; Johnston et al., 2011; Mankin 
et al., 1999; Michael Baker Jr. Inc, Aqua Terra Consultants, & Dynamic 
Solutions LLC, 2015; Privette et al., 2015; Shabani et al., 2017; Sutula et al., 
2016; Yue & Derichsweiler, 2005). 


4.3 Integrated watershed and groundwater 
quality modelling 


To evaluate the impacts of climate and land-use changes on water resour- 
ces (surface and groundwater; quantity and quality) at a regional to water- 
shed scale requires an integration of watershed, groundwater, and receiv- 
ing water quality models which is capable of simulating all the important 
processes of hydrogeological cycle. Table 11 outlines several studies that 
attempted to integrate a watershed model with groundwater and receiving 
water quality models through a sequential coupling. Only few systems are 
reported with the dynamic coupling of all groundwater, watershed, and 
complex receiving water and transport models, including MIKE SHE/ 
MIKE 11 model, which has been described in previous sections. 

Klammler et al. (2013) implemented a sequential coupling of the 
one-dimensional unsaturated water flow and nitrate transport model 
SIMWASER/STOTRASIM with the two-dimensional saturated approach 
of FEFLOW to simulate the nitrate leaching from the soil zone into the 
aquifer Westliches Leibnitzer Feld in Austria to evaluate the impact of 
agricultural practices on groundwater quality. The results of the unsatur- 
ated water model (water and nitrate flux) are provided as the upper time 
series boundary condition to the FEFLOW model. 

Narula and Gosain (2013) applied SWAT, MODFLOW, and MT3DMS 
to model hydrology, groundwater recharge, and non-point nitrate load- 
ings in the Himalayan Upper Yamuna basin in India. The groundwater 
recharge and nitrate (NO,) loads simulated by the SWAT model are linked 
to the groundwater flow model (MODFLOW) and the multi-species 
transport model (MT3DMS). The hydrologic terms simulated by SWAT 
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for each sub-basin were transformed to the system of units specified for 
MODFLOW’s simulation. Grid cells of MODFLOW were associated with 
the geographical extent of sub-basins simulated by SWAT. Groundwater 
limits for the model correspond to those of the surface water basin. These 
boundaries were designated as no flow boundaries. A similar integrated 
modelling framework was developed by Pulido-Velazquez et al. (2015) for 
the integrated assessment of the impact of climate and land use changes on 
groundwater quantity and quality in the Mancha Oriental system in Spain 
by sequentially coupling a watershed agriculturally based hydrological 
model (SWAT) with a groundwater flow model developed in MODFLOW, 
and with a nitrate mass-transport model in MT3DMS. SWAT model out- 
puts (mainly groundwater recharge and pumping, considering new irrig- 
ation needs under changing evapotranspiration and precipitation) were 
used as MODFLOW inputs to simulate changes in groundwater flow, stor- 
age, and impacts on stream-aquifer interactions. SWAT and MODFLOW 
outputs (the nitrate load from SWAT and groundwater velocity field from 
MODFLOW) are used as MT3DMS inputs for assessing the fate and trans- 
port of nitrate leached from the topsoil. 

Ameli and Creed (2017) developed a linked subsurface-surface mod- 
el to assess the continuum of time and distance variations of hydrologic 
connectivity of wetlands in the Beaverhill watershed, Alberta, character- 
ized by the high density of geographically isolated wetlands. A three-di- 
mensional steady-state groundwater-surface water interaction model was 
used to simulate watershed-scale subsurface flow and velocity fields as well 
as to calibrate the infiltration rate. These model results were then used 
to map watershed-scale subsurface connections. The two-dimensional 
transient fill-and-spill surface flow routing approach within the numer- 
ical, physically based HydroGeoSphere model was linked to the output of 
the groundwater model and used to simulate the watershed-scale surface 
water level and overland flow routing, and ultimately to determine the 
surface connectivity of wetlands using a transient water particle tracking 
scheme. The performance of the model was also assessed using chemical 
(Ca, Mg, EC and TDS) and isotopic ('8O and °H) tracer data. 
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Table 11. Examples of Integrated Watershed and Groundwater Quality 


Modelling Studies 


Coupling Author Application Watershed Groundwater Coupling 
region/ model model tool & data 
watershed exchanges 

Sequential Klammler Westliches SIMWASER/ FEFLOW A specific 

et al., 2013  Leibnitzer STOTRASIM add-in module 
Feld, for FEFLOW is 
Austria, 44 developed to 
km? link recharge 
and nitrate 
concentration. 
Narula & Himalayan SWAT MODEFLOW User-specified 
Gosain, Upper + MT3DMS transform 
2013 Yamuna and coupling 
basin, India, for recharge 
11,600 km? and nitrate 
concentration. 
Pulido- Mancha SWAT MODEflow + User-specified 
Velazquez Oriental MT3DMS transform 
et al., 2015 system, and coupling 
Spain for recharge, 
pumping flow, 
and nitrate 
concentration. 
Ameli & Beaverhill A 3D HGS User-specified 
Creed, watershed, ground- transform 
2017 Alberta water-sur- and coupling 
face water for recharge, 
interaction pumping flow, 
model and chemical 


and isotopic 
tracer concen- 
trations 


4.4 Integrated groundwater and receiving water 
quality modelling 


There is a transition zone where groundwater and surface water interact. 


This is an ecologically active zone where contaminants from upland areas 
that are transported by groundwater can be retained within sediments 


and transported to the receiving surface water. Similarly, contaminants 


discharged to the surface water can be a source of contamination to 
groundwater if the surface water recharges the underlying aquifer. Both 


sediments and surface water provide a pathway by which contaminants 
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can enter to the groundwater systems. The transition zone is strongly 
influenced by the dynamic exchanges between groundwater and surface 
water and changing biogeochemical conditions (Bobba, 2012). 

Traditional approaches to model groundwater-surface water inter- 
actions often focus on representing one hydrological system in detail and 
the other as a boundary condition without explicitly considering the ef- 
fects of feedback between the two systems. Models in this category are 
integrated using a sequential coupling approach. For example, to simulate 
pit lake water quality, modelling knowledge from different scientific do- 
mains such as groundwater, lake circulation, hydrochemistry, and lim- 
nology needs to be combined. The modelling system MODGLUE couples 
the groundwater flow and transport model PCGEOFIM with the lake 
circulation water quality model CE-QUAL-W2 and the hydrochemical 
model PHREEQC (Müller et al., 2008). Jia et al. (2015) developed a linked 
surface water and groundwater simulation model to assess the impact of 
a trans-basin water diversion project on the groundwater. By using results 
of the surface water simulation as input for the groundwater simulation, a 
surface water quality WASP and a groundwater model MODFLOW plus 
MT3 were sequentially coupled to simulate the water levels and four con- 
taminants (NH3-N, COD, Mn, F, As). 

In cases where relatively large, dynamic, bidirectional exchanges are 
anticipated, the interfacial processes cannot be adequately represented by 
using an uncoupled interaction to represent surface water in a groundwater 
model or with a source term to represent groundwater flux in a surface 
water model. Although the existing ready-to-use integrated subsurface and 
surface water and solute transport models such as HGS and MIKE SHE 
are capable of simulating groundwater-—surface water interactions, they do 
not appear to provide a full representation of sediment/benthic processes. 
Mugunthan et al. (2017) developed an interface module that holistically 
simulates fate and transport by dynamically coupling two commonly used 
models, AQFATE and SEAWAT, to simulate surface water and ground- 
water hydrodynamics, while providing an enhanced representation of the 
processes in the transition zone. AQFATE is an enhanced version of the 
EFDC model (Connolly et al., 2000). SEAWAT is a groundwater model 
developed by USGS which combines MT3DMS'’s solute transport capabil- 
ities with MODFLOW to simulate density effects on groundwater flow 
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(Langevin et al., 2008). The interface code is developed in FORTRAN, the 
same language used in the original SEAWAT and AQFATE models. At each 
groundwater sub-model timestep, the modelling framework represents the 
surface water body as a boundary condition. Constituent concentrations 
(temperature, salinity, or contaminants) are passed to SEAWAT. Upon 
completion of the first and subsequent groundwater sub-model timesteps, 
the AQFATE sub-model is simulated for the corresponding period with the 
flows and mass fluxes calculated by SEAWAT at the interfacial grid cells 
passed to AQFATE through intermediate variables in the interface module 
code. The modelling framework was tested with a published test problem 
and applied to evaluate field-scale two- and three-dimensional contam- 
inant transport. The model accurately simulated concentrations of salin- 
ity from a published test case. Table 12 lists some examples of integrated 
groundwater and receiving modelling studies. 


Table 12. Examples of Integrated Groundwater and Receiving Water 
Quality Modelling Studies 


Coupling Author Application Groundwater Receiving Coupling 
region/ model water tool & data 
watershed quality exchanges 

model 

Dynamic Mugunthan A former oil SEAWAT AQFATE (an | User-devel- 

et al., 2017 refinery site enhanced oped inter- 
in Western version of face module, 
Canada EFDC) coupling 
flow, and 
concentra- 
tions 
Sequential Muller etal., Severalmine PCGEOFIM CE- User- 
2008 pit lakes in QUAL-W2 developed 
Germany PHREEQC interface, 
coupling 
flow, water 
level, and 
water quality 
parameters 
H. Jia et al., Chaobai MODFLOW + WASP User defined 
2015 River alluvial MT3D text files, 
plain, Beijing, coupling 
China flow, water 
level, con- 
centrations 
of NH3-N, 
COD, F, As 
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4.5 Integrated atmospheric deposition and water 
quality modelling 


Atmospheric wet and dry deposition can be important non-point-source 
contributors to total pollutant loadings to water bodies, both through dir- 
ect deposition to water bodies and deposition to watersheds with subse- 
quent transport into water bodies. In a study of the nitrogen budgets of 
16 catchments in the northeastern United States, atmospheric deposition 
was found to be the largest source of nitrogen input to the catchments, 
contributing about 31% to the overall budget (Boyer, 2002). Atmospheric 
deposition can affect ecosystems in numerous ways including acidifica- 
tion and eutrophication. Acidification of lakes and streams is primarily 
caused by the atmospheric deposition of sulfur (S) and reactive nitrogen 
(N) to watersheds, with some impact from direct deposition to lakes. The 
deposited chemicals undergo subsequent biogeochemical cycling and the 
transfer of chemicals to surface water systems (Paerl, Dennis, & Whitall, 
2002; Sullivan et al., 2008). 

Quantification of the atmospheric deposition is important to water 
quality studies. Watershed-scale fate and transport models such as SWAT 
and HSPF use this information to estimate loadings to rivers and water- 
sheds, for use in TMDL developments and other water quality assessment 
and management plans. However, obtaining good estimates of atmos- 
pheric wet and dry depositions can be challenging. Direct measurement 
of deposition, particularly dry deposition, can be difficult and very expen- 
sive to monitor at several sites in a watershed (Schwede, Dennis, & Bitz, 
2009). Atmospheric deposition models, generally classified as Eulerian 
and Lagrangian models, can be used to fill in spatial or temporal holes left 
by a monitoring program and predict future conditions due to growth or 
regulatory changes (NEIWPCC 2017). Eulerian models perform calcula- 
tions of atmospheric chemistry, transport, and deposition of pollutants 
based on grids. Eulerian models are effective for capturing the complex 
nonlinear chemistry necessary to model ozone, nitrogen, sulfur, and 
mercury accurately. Examples of Eulerian models include the Regional 
Acid Deposition Model (RADM), the Regulatory Modelling System for 
Aerosols and Deposition (REMSAD), and the Community Multiscale 
Air Quality (CMAQ) model. Lagrangian models generally work well for 
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toxic compounds that have simple decay or linear atmospheric chemistry. 
These models track emission plumes that spread out toward some recep- 
tors, such as an estuary, where deposition is taking place, based on the 
receptor’s chemical and physical parameters and meteorology. Examples 
of Lagrangian models include the Regional Lagrangian Model of Air 
Pollution (RELMAP) and the California Puff Model (CALPUFF). 

For integrated air and water quality management, there have been 
significant advances in the development of integrated airshed, watershed, 
and water body modelling and analysis technologies. A sequential coup- 
ling approach is generally used to link the output of a deposition model to 
watershed and receiving water quality models. In the United States, mod- 
els of the Chesapeake Bay airshed, watershed, and tidal waters have been 
created and linked to model daily atmospheric deposition loading and 
the impacts on bay water quality and resources (e.g., underwater grasses, 
benthic communities, pelagic fish habitats) (Ackermann, 1997). In par- 
ticular, the Regional Acid Deposition Model (RADM) has been used to 
delineate the airshed contributing nitrate to the Chesapeake Bay water- 
shed and water surface (Dennis, 1997). Burian et al. (Burian et al., 2002) 
developed an integrated modelling framework composed of a CIT urban 
air chemistry model, a SWMM urban runoff model, and a WASP water 
quality model. The models were linked to simulate the fate and transport 
of air emissions of nitrogen compounds in the air, urban watersheds, sur- 
face water runoff, and a coastal receiving water body. The model linkage 
is demonstrated by evaluating the potential water quality implications of 
reducing NOx emissions by 32%, volatile organic compound emissions by 
51%, and ammonia emissions by 30%, representing changes from the 1987 
levels to the proposed 2000 target levels in Los Angeles, California. 

Sequential coupling requires post-processing the dry and wet depos- 
ition outputs from a deposition model into an input format required by 
watershed and receiving water quality models. For example, a tool called 
the Watershed Deposition Tool was developed for providing the linkage 
between air and water quality modelling and for analyzing related non- 
point-source impacts on the watershed. Using a gridded output of atmos- 
pheric deposition from the CMAQ model, the tool calculates the average 
per unit area and total deposition to selected watersheds and sub-water- 
sheds (Schwede et al., 2009). 
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4.6 Integrated load allocation and water quality 
modelling 


Load allocation is the distribution of pollutant loadings among point and 
non-point sources in a watershed such that the receiving water body is en- 
sured to be compliant with water quality standards. In the United States, 
Section 303(d) of the Clean Water Act established the total maximum 
daily load (TMDL) approach to water quality management. A TMDL is 
the maximum loading rate of a pollutant that can be sustained in a water 
body without water quality impairments. It also specifies the maximum 
amount of a pollutant that a waterbody can receive and still meet water 
quality standards and allocates pollutant loadings among point and non- 
point pollutant sources. A TMDL is the sum of the individual waste-load 
allocations (WLAs) for point sources, load allocations (LAs) for non-point 
sources, and the natural background with a margin of safety (USEPA, 
2008). 

In most load allocation analyses or TMDL developments, load reduc- 
tion scenarios are evaluated through trial and error (TAE) simulations, 
using a process-based watershed water quality model. However, the TAE 
simulation method does not necessarily generate cost-effective, reliable, 
and equitable load allocations (Y. Jia & Culver, 2006). To overcome the 
limitation of TAE scenario analysis, coupled simulation-optimization 
models have been built and recommended to develop optimal load alloca- 
tions (USEPA, 2008). One of the challenges of applying optimal manage- 
ment strategies for load allocation is defining an appropriate objective 
function, decision variables, and constraints. Potential objectives may be 
maximization of equity among sources, minimization of total load reduc- 
tion, maximization of total net benefit, or minimization of total cost. For 
example, the same percentage reduction in load contribution could be ap- 
plied to sources of similar types, or the maximum difference in percentage 
reduction among sources can be defined as constraints. Although it may 
be desirable to maximize the benefit or minimize the cost of load reduc- 
tions, this requires substantial site-specific information about the types 
of management alternatives and associated costs and benefits. Multiple 
objective optimization (Allam et al., 2016) and game theoretical models 
(Nikoo, Beiglou, & Mahjouri, 2016) have often been developed for load 
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allocation in order to take into account stakeholders in the negotiation 
processes. 

The water quality constraints are generally based on simulated receiv- 
ing water quality and water quality guidelines or limits. The optimiza- 
tion-simulation load allocation model is essentially an iteratively coupled 
model, which requires repeated calls of the water quality model to solve 
the optimization model. For instance, multiple objective models were de- 
veloped and applied to the Gharbia drain in the Nile Delta in Egypt during 
both the summer and winter seasons of 2012 through the integration of 
a water quality model QAUL2Kw and a genetic algorithm, by considering 
the total waste load abatement and the inequity among waste dischargers. 
The steady water quality model was directly coupled and iteratively run 
to produce BOD, concentrations for a particular arbitrary treatment level 
during the process of optimization. 

To reduce the prohibitive computational resources often required by 
the iteratively coupled model, an indirect simulation-optimization frame- 
work utilizing a response matrix approach is commonly used to replace 
the iterative calls of the water quality model (Y. Jia & Culver, 2006). In 
this approach, output of the water quality model runs are used to derive a 
linear response matrix (Y. Jia & Culver, 2006), load delivery factors (Shortle 
et al., 2016), transfer coefficients (Nikoo et al., 2016; Zolfagharipoor & 
Ahmadi, 2016), or nonlinear stressor-response relationships (F. Zhou et 
al., 2015), which are then read by the optimization model to find the solu- 
tion. Therefore, an iterative coupled optimization-simulation problem is 
transformed into a sequentially coupled simulation-optimization model. 
A robust optimization approach to minimize total load reduction was 
successfully developed by Jia and Culver (2006) using sequential coupling 
with a response matrix and applied to the fecal coliform TMDL study in 
the Moore’s Creek watershed located in Albemarle County, Virginia, USA. 
The response matrix was defined with elements representing temporal 
changes in water quality with unit load reduction for each source derived 
based on water quality simulation time series results. One advantage of 
the robust formulation of TMDL allocations is that the uncertainty of the 
watershed simulation model, HSPF, is incorporated into the load alloca- 
tion optimization model by introducing the probability of acceptable par- 
ameter sets of the watershed model and corresponding simulated baseline 
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Table 13. Examples of Integrated Load Allocation and Water Quality 


Modelling Studies 


Coupling Author Application Load Water Coupling 
region/ allocation quality tool & data 
watershed model model exchanges 

Iterative Allam et al., Gharbia Minimize QUAL2Kw Direct 

2016 catchment, total waste coupling 
Egypt, 2940 abatement; steady water 
km? minimize quality model 

inequity (BOD5/ 
among DO) with 
wastewater optimization 
dischargers model 

Sequential Y. Jia & Moore’s Creek Minimize to- HSPF Response 

Culver, watershed, VA tal weighed matrix for 

2006 load reduc- instream fecal 

tion coliform con- 
centrations to 
1% reduction 
of loads 

Nikoo et Zarjub River, Non- QUAL2Kw Transfer 

al., 2016; Iran cooperative coefficients 

Zolfaghar- and and trading 

ipoor & cooperative ratios, 

Ahmadi, game determined 

2016 theoretic based on the 

multiple- results of a 

pollutant calibrated 

waste load QUAL2Kw 

allocation model for 

models BOD/DO and 
TN 

Zhou et al., Swift Creek Enhanced-in- CE- Nonlinear 

2015 Reservoir, terval linear QUAL-W2 stressor- 
Chesterfield program- response 
County, VA ming for nu- relationships 

trient TMDL 
allocation 

Shortle et Chesapeake Cost minimi- Chesa- Delivery 

al., 2016 Bay watershed, zation static peake Bay factors, land 
MD and dynamic Watershed areas, and 

optimal Model baseline 
models (CBWM) nutrient 
loadings 
based on 
watershed 
modelling 
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source load contributions into the objective function and water quality 
constraints, respectively. A total of 381 acceptable parameter sets for the 
Moore’s Creek HSPF model were established using the Monte-Carlo- 
based generalized likelihood uncertainty estimation (GLUE) approach 
with 50,000 HSPF runs (Beven & Binley, 1992; Y. Jia & Culver, 2008). The 
likelihood value of each of these is calculated using a fuzzy logic proced- 
ure. The robust optimization model was then solved using a genetic algo- 
rithm. Some of the integrated load allocation and water quality modelling 
studies are listed in Table 13. 


4.7 Integrated water allocation and water quality 
modelling 


Water allocation is the combination of actions that enable water users to 
take or to receive water for beneficial purposes according to a recognized 
system of rights and priorities (UN-ESCAP 2000). Water allocation is cen- 
tral to the management of water resources, which often engages multiple 
stakeholders with conflicting interests. Dinar et al. (1997) discuss four 
basic institutional mechanisms for water allocation: user-based allocation, 
marginal cost pricing, public allocation, and water markets allocation. 
Water allocation models have been developed for different purposes such 
as water rights allocation (Labadie, 1995; L. Wang, Fang, & Hipel, 2007), 
economic optimal water allocation (McKinney, 1999), and cooperative, 
fair, efficient, and sustainable water allocation (L. Wang, Fang, & Hipel, 
2008). Although the inseparable interaction of water quantity and quality 
clearly exists in all river basins, most water allocation models focus on 
water quantity with interactions, if any, accounted for by superficial trial 
and error processes. This trial and error water allocation with considera- 
tion of water quality is achieved based on a simple sequential coupling, i.e., 
a water allocation model is run to provide flow inputs for subsequent water 
quality modelling (Salla et al., 2014). To eliminate the limitations of the 
trial and error approach, water quality models have been directly linked to 
optimal water allocation models and are iteratively called to run for each 
potential water allocation. One typical iteratively coupled water allocation 
model is the MODSIMQ model developed by Dai and Labadie (2001). The 
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Table 14. Examples of Integrated Water Allocation and Water Quality Modelling Studies 


Coupling Author Application region/ Water Water Coupling tool & data 
watershed Allocation model quality model exchanges 
Iterative Dai & Labadie, Arkansas River MODSIMQ QUAL2E Frank-Wolfe nonlinear 
2001 Basin, CO programming algorithm, cou- 
pling salinity concentrations 
D. Liu et al., Northwest Pearl Multi-objective Water 1D advection- Non-dominated sorting GA-II 
2013 River Delta, China quantity and waste dispersion water (NSGA-II) algorithm, coupling 
load allocation model: quality model COD concentration 
minimize water shortages, 
maximize economic 
interest, maximize waste 
load discharges subject to 
water quality targets 
Sequential Salla et al., Araguari River SIMGES module of GESCAL module Text files, coupling oxygen, 
2014 basin, Brazil AQUATOOL of AQUATOOL biochemical oxygen demand, 


organic nitrogen, ammonia, 
nitrate, and total phosphorus 


Tavakoli et al., 
2014 


Dez River, Iran 


An optimization model 
based on equitable 
allocation of water to 
users in proportion to their 
water demands 


Soil,Water, Atmo- 
sphere, and Plant 
(SWAP) simulation 
model 


Iterative linear programming 
CILP), coupling 5 meta-mod- 
els, each zone representing 
the relationships between 
quantity and quality (TDS) of 
return flow versus the allocat- 
ed water 


Heydari et al., 
2016 


Zayanderood river 
basin, Iran 


Multi-Objective 
Optimization Model: 
minimize relative water 
deficit, minimize annual 
groundwater level 
changes, minimize the 
groundwater quality 
change 


MODFLOW and 
MT3DMS models 


Two surrogate models, name- 
ly an Artificial Neural Net- 
work model for groundwater 
level simulation and a Genetic 
Programming model for TDS 
concentration prediction 
were coupled with NSGA-II 


MODSIM river basin water rights planning model developed at Colorado 
State University is extended by MODSIMQ, integrating with the Frank- 
Wolfe nonlinear programming algorithm to directly include conservative 
routing of water quality constituents, maintenance of salinity load mass 
balance, and the imposition of constraints on water quality concentra- 
tions. Water quality constraints can be imposed based on (i) quality stan- 
dards for certain river reaches, (ii) irrigation water quality control, (iii) 
water quality preference for demand nodes, and (iv) groundwater quality 
rehabilitation. An iterative procedure between MODSIMQ and QUAL2E 
assures convergence to solutions that satisfy water right priorities, while 
attempting to maintain minimum streamflow and water quality require- 
ments. In the literature, there are a few studies on simultaneous water re- 
sources and waste load allocation in river basins, in which waste loads are 
also included as decision variables and objective functions (D. Liu et al., 
2013). 

Generally, integration of a nonlinear simulation model in a manage- 
ment model is difficult and computational time to achieve the optimal 
solution may be a constraint (Singh, 2014). The required computational 
time can be reduced via approximations of the simulation model by using 
simplified response matrixes or surrogate models as alternatives to actual 
complex numerical models (Heydari, Saghafian, & Delavar, 2016; Tavakoli 
et al., 2014). The actual complex water quality models are sequentially 
coupled to water allocation models through the approximate relation- 
ship or surrogate models between water quality and quality. Table 14 lists 
some examples of integrated water allocation and water quality modelling 
studies. 
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5.0 MODELLING IN THE ATHABASCA 
RIVER BASIN - CASE STUDY 


The Athabasca River, located in Alberta, Canada, originates at the 
Columbia Ice Fields near the Alberta—British Columbia border and flows 
approximately 1300 km northeast before entering Lake Athabasca at the 
northeastern corner of Alberta (Figure 8). Water from Lake Athabasca 
flows into the Slave River and joins the Mackenzie River, which eventually 
enters the Arctic Ocean. The elevation of the watershed varies from more 
than 3000 m a.s.l. in headwaters in the Columbia Icefield to about 205 m 
a.s.l. at its outlet in Lake Athabasca. The Athabasca River basin is physic- 
ally and ecologically diverse and covers an area of approximately 159,000 
square kilometers. The region is endowed with many natural resources 
such as forests, coal, minerals, agriculture, and oil and gas. The Athabasca 
oil sands are large deposits of bitumen or extremely heavy crude oil and 
are the largest reservoir of crude bitumen in the world and the largest of 
three major oil sands deposits in Alberta (along with the nearby Peace 
River and Cold Lake deposits). 

Interest in bitumen from the oil sands dates back to the 1930s, when 
intensive industrial development of the region began, and reached a peak 
in the late 1960s. In the mid 1990s, new project applications and the expan- 
sion of existing oil sands operations were underway. As a result, the lower 
Athabasca watershed was getting more attention for CEA in response to 
the large-scale oil sands operations that can disturb environmental condi- 
tions in the watershed. In-situ and open pit mining developments can dir- 
ectly affect subsurface and surface hydrological processes such as runoff, 
soil moisture, and infiltration. For instance, the utilized steam-assisted 
gravity drainage (SAGD) has a potential to modify the hydrogeological 
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regime of the basin. A part of the surface water used in SAGD operations 
(i.e., steam production and then injection) is lost to the bitumen recov- 
ery. (It occupies the space previously occupied by oil in formations.) In 
response to the changes in environmental conditions, the Alberta govern- 
ment took steps to initiate several strategies and plans. In the late 1990s, 
the Alberta government designed its Regional Sustainable Development 
Strategy (RSDS) to address potential cumulative environmental effects 
in the Athabasca oil sands area. The aim of the RSDS was to provide a 
framework for managing cumulative environmental effects and to ensure 
sustainable development in the Athabasca oil sands area. The strategy pri- 
oritized 72 environmental issues which had been divided into a list of 14 
themes and 3 priority categories that should be assessed in response to the 
oil sands development. 

In 2008, the Government of Alberta commenced a comprehensive 
initiative to develop a new approach for managing cumulative effects by 
releasing a land-use framework, known as the Lower Athabasca Regional 
Plan (LARP). The LARP is a comprehensive, forward-thinking, and legal- 
ly binding roadmap that enhances environmental management, addresses 
growth pressures, and supports both human and ecosystem needs, while 
balancing social, environmental, and economic outcomes. The regional 
plan considers the cumulative effects of all activities on air, water, and bio- 
diversity. Five environmental management frameworks were developed 
under the Lower Athabasca Regional Plan, including the Air Quality, 
Surface Water Quality, Groundwater, Surface Water Quantity, and 
Tailings Management frameworks. A biodiversity management frame- 
work is being developed. 

In 2012, the Governments of Alberta and Canada embarked on a 
new plan, known as the Joint Oil Sands Monitoring Program (JOSM), 
to ensure a comprehensive and systematic monitoring and reporting of 
environmental conditions in the lower Athabasca River basin to support 
sustainable resource development. 
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Figure 8. Athabasca River Basin (ARB) 
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The above-mentioned programs, plans, and strategies were composed 
of different themes and sub-themes to build an understanding of the en- 
vironment as a whole and to ensure that any major environmental impacts 
are considered. As a result, many individual studies were conducted with 
a focus being on either theme or environmental processes. Some of the 
environmental modelling studies in the region (especially in the oil sands 
areas) are described in the following sections. 


5.1 Hydrodynamic and water quality modelling in 
the Athabasca River 


Water quality modelling of the Athabasca River began in 1984 when the 
first model was implemented using Water Quality for River and Reservoir 
Systems (WQRRS) for the main stem from Hinton to Lake Athabasca 
(Charles Howard and Associates Ltd., 1984). The calibration process was 
challenging due to limited monitoring data (e.g., hydraulic, non-point and 
point source loadings, and in-stream water quality data). Flows under 
ice cover could not be simulated. A dissolved oxygen simulation model 
named DOSTOC (Dissolved Oxygen STOChastic model) had been widely 
used since its development in 1987, originally for the Planning Division of 
Alberta Environment (HydroQual Consultants Inc. & Gore and Storrie 
Ltd., 1989). In 1988, trial runs of the DOSTOC model were undertaken for 
the Athabasca River using data collected during the 1987 and 1988 win- 
ters (Culp & Chambers, 1994). The model was used to simulate the water 
quality in the Athabasca River under various levels of development in the 
basin and low flow conditions during ice-coved periods. The calibration 
of the DOSTOC model was updated with data collected during winter 
synoptic surveys in 1988, 1989 (Macdonald & Hamilton, 1989), and 1990 
(Macdonald & Radermacher, 1992), and then further evaluated and valid- 
ated by the Northern River Basins Study (NRBS) using winter survey data 
obtained in 1991 and 1992 (Macdonald & Radermacher, 1993) as well as 
1993 and 1994 (Chambers et al., 1996; Pietroniro, Chambers, & Ferguson, 
1998). The DOSTOC model is a component of the Stochastic River Quality 
Model (SRQM) that consists of three modules: DOSTOC for dissolved oxy- 
gen, NUSTOC for nutrients, and UNSTOC for user-specified substances. 


120 5.0 Modelling in the Athabasca River Watershed - Case Study 


The model was formulated as a one-dimensional steady state model and 
based on the analytical solution to the stochastic version of the Streeter- 
Phelps equation, assuming that random parameters follow normal distri- 
butions. The model can be run in either deterministic or stochastic mode. 
Hydraulic characteristics of the river were represented by the Leopold- 
Maddock equations, which are exponential functions relating mean 
depth, top width, and mean velocity to river discharge. Ice cover processes 
were not simulated, and effects were implemented by assuming zero vola- 
tilization and reduced photolysis and biodegradation rates. Water quality 
processes represented in the model include atmospheric reaeration, decay 
of BOD, and nitrogenous oxygen demand (NOD) in the water column, 
photosynthesis and respiration, and benthic SOD (McCauley, 1997).The 
calibrated model was used to evaluate the impacts of pulp and paper in- 
dustry development in the region. It was determined that although the 
treated pulp effluents from two mills (Wildwood at Hinton and Millar 
Western at Whitecourt) would increase BOD and ammonia concentra- 
tions directly downstream of the mills, the impact is negligible further 
downstream on the river. As the DOSTOC model is a steady state model 
over a short period, it could not capture the cumulative effects over time 
from spatially distributed sources. 

During the NRBS, one-dimensional dynamic models with separate 
and interacting water column and bed sediment compartments were also 
developed to simulate the fate and transport of organic chemicals for 
the Athabasca (Hinton to Old Fort) and Wapiti/Smoky Rivers, using the 
USEPA’s Water Quality Analysis Simulation Program version 4 (WASP4) 
(Golder Associates Ltd., 1997a, 1997b). The Leopold-Maddox method was 
added to WASP4 code so that water column velocities, cell volumes, and 
mass exchange areas would be updated at each timestep for the water col- 
umn segments. This approximate approach is suitable for gradually varied 
flows. A sediment transport algorithm for the Athabasca River, developed 
by Krishnappan et al. (1995) was incorporated in WASP4 to predict resus- 
pension and deposition rates within the Athabasca River contaminant fate 
model. Out of the seven selected organic chemicals, the best calibration 
over the two-year period of 1992-1993 was achieved for 2,3,7,8-TCDF. For 
other substances, especially phenanthrene, observed data were so sparse, 
or conflicting, that it was not possible to evaluate the calibration. By 
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incorporating the sediment transport algorithm into WASP4, the model 
predicted a very dynamic exchange of sediment between the water col- 
umn and bed, with an accumulation of fine bed sediment over the late fall 
and winter, removal by resuspension during the spring freshet, and very 
little net accumulation during the summer. 

Subsequently, a one-dimensional model, WASP7 (water quality 
analysis simulation program), was set up using a kinematic wave rout- 
ing scheme and was hydrodynamically calibrated and validated for the 
Lower Athabasca River with 1999-2008 data (Kannel & Gan, 2013). The 
model represented the field data quite well except during winter seasons 
(mid-November to April). The model was applied to investigate the poten- 
tial impact of oil sands processed water (OSPW) in the event that OSPW, 
which contained naphthenic acids (NAs), was accidentally discharged to a 
stretch of the Athabasca River, simulating NAs as a lumped state variable 
and assuming it is degraded by natural dilution, biodegradation, sorption, 
photodegradation, or combinations of these processes. NAs in the Lower 
Athabasca River were predicted to be sensitive to changes in the discharge 
rate and concentrations of OSPW, as well as the rates of photodegradation 
and biodegradation. 

Numeric modelling of flow and transport processes in the Athabasca 
Oil Sands Region is challenging due to the complex morphology, cold 
climate, and highly variable flows in the river. Being located in northern 
Alberta, the Lower Athabasca River has some special hydrodynamic and 
water quality characteristics, especially the ice formation, jams, melting 
and break-up processes, and the relatively long ice-cover period for the 
river. Early attempts to model flow in the Lower Athabasca River used 
one-dimensional models that were implemented based on approximate 
and simplified rectangular cross-sections to represent channel geometry. 
Khanna and Herrera (Khanna & Herrera, 2002) applied the cdg1-D model 
in the Lower Athabasca River basin to estimate high flows during open- 
water season. The cdg1-D model, originally developed at the University of 
Alberta, solves the St. Venant equations by the finite element method using 
the characteristic dissipative Galerkin (cdg) scheme. Trillium Engineering 
and Hydrographics Inc. (Trillium Engineering and Hydrographics Inc., 
2003) conducted field surveys to measure transverse mixing coeffi- 
cients and travel time in the Lower Athabasca River during ice-covered 


122 5.0 Modelling in the Athabasca River Watershed - Case Study 


winter low flow periods at five key locations (short reaches) downstream 
of the Firebag River, Ells River, Muskeg River, Steepbank River, and Fort 
McMurray. A HEC-RAS model for each reach was constructed to evaluate 
the hydraulic roughness and to predict water levels and mean velocities for 
a range of discharges. Two-dimensional flow characteristics required for 
the mixing analysis were evaluated using a lateral discharge distribution 
approach. Two methods of evaluating transverse mixing coefficients were 
employed: (i) an analytical model with reach-average hydraulic charac- 
teristics and (ii) a numerical model, TRSMIX, which employed local hy- 
draulic characteristics. These dimensionless coefficients may be applied 
over the range of typical winter discharges that occur as long as the river 
is ice-covered. 

Hydrodynamic and fish habitat modelling have been performed for 
reaches of the Lower Athabasca River and Peace-Athabasca Delta using 
the River2D model calibrated with synoptic bathymetry and hydrometric 
survey data obtained during summer and winter (ice-covered) periods as 
part of a multi-year program led by the Surface Water Technical Group 
of the Cumulative Environmental Management Association (CEMA) 
(AMEC-nhc, 2009). The program aimed to assess in-stream flow needs 
and evaluate fish habitats for open water and winter conditions for the 
Lower Athabasca River, in which five study segments along the Lower 
Athabasca River were investigated, including: Reach #1 - Athabasca Delta 
below Embarras, Reach #2 - Embarras, Reach #3 - Poplar Point, Reach 
#4 - Bitumount, and Reach #5 - Northlands. River2D is a two-dimen- 
sional finite element-based numerical hydrodynamic model that was de- 
veloped at the University of Alberta to simulate the depth-averaged flow 
characteristics in a river segment (Steffler & Blackburn, 2001). To match 
measured velocities and water levels, a range of roughness heights was 
adopted for the calibration of flows: 3.0 mm for sand, 500 mm for the 
cobble regions, and 10 mm for ice at Reach #4 (Trillium Engineering and 
Hydrographics Inc., 2004, 2005); 10 mm for sand and 150 mm for ice at 
Reach # 2 (Northwest Hydraulics Consultants Ltd., 2007a), and 1 mm for 
sand and 150 mm for ice at Reach #3 (Northwest Hydraulics Consultants 
Ltd., 2007b). Katopodis and Ghamry (2005) conducted similar ice-cov- 
ered hydrodynamic model calibrations and comparisons for three reaches 
of the Lower Athabasca River (Fort McKay below Peter Lougheed Bridge, 
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Bitumount, and Northlands). It was found that the applied bed and total 
roughness along the thalweg profiles for Bitumount and Northlands 
Reaches were comparable for similar substrate sizes. The Northlands 
Reach has the coarsest bed materials over most of the thalweg profile, 
the Peter Lougheed Reach has the finest ones, and the Bitumount Reach 
is in-between. Although the applied ice roughness differed between the 
three reaches, its low values had a small effect on the composite roughness. 

In the Peace-Athabasca Delta, 2D-River models were developed and 
calibrated for the two divergence areas under each of the open water and 
under-ice conditions. The ice thickness data obtained in the winter survey 
was processed to produce River2D ice input files (AMEC-nhc, 2009). A 
Riverl1D model for the Peace-Athabasca Delta was also developed at the 
University of Alberta (Andrishak & Hicks, 2009, 2011), for the primary 
purpose of simulating river discharges across the Lower Athabasca Region 
to provide boundary conditions for more detailed (e.g., River2D) hydraul- 
ic and habitat modelling at flow divergence sites. Both 1D- and 2D-River 
models for the Peace-Athabasca Delta were updated in 2014 with new 
survey data (Hatfield Consultants, 2014). Note that the 1D- and 2D-River 
models for the Lower Athabasca River and Peace-Athabasca Delta were 
all calibrated with steady-state runs under steady-state conditions and 
specified temporally constant ice thickness inputs interpolated from data 
collected from synoptic surveys. Steady-state runs were performed for fish 
habitat modelling and assessment due to their limited functionality for 
unsteady modelling. Efforts were made and preliminary results were ob- 
tained by updating the 1D-River (Andrishak et al., 2008) and 2D-River 
model (Wojtowicz et al., 2009) to include thermal ice processes to simulate 
the freeze/thawing of the Athabasca River. 

Integrated hydrodynamic and water quality modelling over a con- 
tinuous period covering open-water and under-ice conditions in the 
Athabasca Oil Sands Region started in the late 2000s. A two-dimensional 
vertically-averaged Environmental Fluid Dynamics Code (EFDC) mod- 
el with hydrodynamic and eutrophication water quality models for the 
Lower Athabasca River from Fort McMurray to Old Fort was developed 
by TetraTech in 2009. The model was used to simulate an eight-year period 
from 2000 to 2007 with a partial calibration due to limited monitoring 
data. The model was enhanced by Dynamic Solutions International LLC 
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(DSI) in a scoping study by updating model bathymetry to the best avail- 
able data and adding preliminary setup of sediment and toxics modules 
without calibration. The model did not simulated the ice formation and 
melting processes but required external inputs for ice cover thickness 
(Dynamic Solutions International LLC., 2012). A two-dimensional lat- 
erally averaged model for the Upper Athabasca River (Hinton to Grand 
Rapids) was developed using CE-QUAL-W2 to simulate the hydrodynam- 
ics, and DO including the ice formation and melting processes over the 
period from 2000 to 2006 (Martin et al., 2013). The modelling results 
indicated that the DO concentration in the Upper Athabasca River was 
very sensitive to the sediment oxygen demand (SOD), which represented 
about 50% of the DO sink in winter. The model was applied under steady- 
state winter low-flow scenarios to predict assimilation capacity for the 
BOD load. The CE-QUAL-W2 model has also been utilized to develop 
for CEMA a general two-dimensional laterally averaged integrated hydro- 
dynamic and water quality model, simulating oil sand pit lakes. Named 
CEMA Oil Sands Pit Lake Model (OSPLM), it can simulate potential water 
quality implications of mature fine tailings placement in pit lakes (Berger 
& Wells, 2014; Golder Associates Ltd. & ERM, 2012; Prakash et al., 2015; 
Vandenberg et al., 2014). 

Modelling the fate and transport of fine sediments and associated 
chemical constituents in the Lower Athabasca River originating from 
natural and potential anthropogenic sources is recognized as a subject 
of increasing importance, as studies have shown that the concentrations 
of sediment-associated chemicals such as PAHs and heavy metals in the 
Lower Athabasca River are affected by development activities (Droppo et 
al., 2018). To quantify and model the sources, transport, and fate of chem- 
icals, a reliable integrated hydrodynamic, sediment transport and water 
quality model of the Lower Athabasca River is needed. Experimental and 
field assessment of sediment dynamics and associated chemicals in the 
Lower Athabasca River and tributaries have been investigated in sever- 
al previous studies under the JOSM program. Droppo and Krishnappan 
(Droppo & Krishnappan, 2016) applied a modelling approach combining 
two existing models (RIVFLOC and MOBED) to simulate the hydro- 
phobic, cohesive sediment transport in the Ells River. Using fine sediment 
transport parameters derived from laboratory flume experiments (e.g., 
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settling velocity of sediment as a function of floc size and the critical shear 
stresses for deposition) and the calculated flow field from the MOBED 
model (using field survey data such as cross-sectional geometry, river 
slope, grain size of bed material, and discharge), the RIVFLOC model was 
used to predict the transport characteristics of the hydrophobic Ells River 
sediments. Although flocculation was shown to occur with increasing floc 
size downstream, there was a breakpoint at approximately 50 um where 
the settling velocity decreased with increasing floc size due to a decreas- 
ing floc density. The high bed shear stresses in the Ells River also negated 
the influence of flocculation on settling. The entrapment process was thus 
concluded as an important aspect of sediment dynamics within high- 
energy cobble/gravel bed rivers, particularly where sediments are hydro- 
phobic like those from the McMurray Formation in Northern Alberta. 

It is well known that integrated watershed and water body model- 
ling systems or tools are needed to support the assessment of cumulative 
effects from climate change, land use changes, developments, and oper- 
ational activities in the oil sands region. Nevertheless, there is no exist- 
ing modelling system or tool that possesses the capability to simulate 
dynamic interactions among environmental and human sub-systems and 
their impacts on water quality and the aquatic habitat health of the Lower 
Athabasca River, accumulated over both spatial and temporal scales. A 
unique, comprehensive, and practical system for integrated watershed 
and water quality modelling in the Athabasca Oil Sands Region was de- 
veloped by Golder Associates, which has been applied to support a series 
of environmental impact assessment (EIA) studies of oil sands develop- 
ment projects (Golder Associates Ltd., 2003a, 2003b, 2004a, 2004b; Teck 
Resources, 2011). The component models include, CALMET/CALPUFF 
air quality dispersion model, regional 3D MODFLOW groundwater 
model, MT3D solute-transport model, regional HSPF hydrologic mod- 
el, CE-QUAL-W2 Pit lake hydrodynamic model, Golder Pit Lake Water 
Quality Model, quasi-dynamic 2D ARM (Athabasca River Model) water 
quality model, steady-state sediment quality model, and habitat suitabil- 
ity models. Component models were calibrated for historical periods. 
Future scenarios and non-point source loadings were estimated based on 
the watershed modelling of LULC changes and development scenarios, 
with the consideration of operational releases from oil sands development 
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as well as withdrawals from the Lower Athabasca River. As HSPF is a 
semi-distributed watershed model, the system does not simulate fully dis- 
tributed interactions between surface and subsurface water. Being a tool 
supporting project EIAs, the system predicts impacts under certain future 
snapshot scenarios and does not simulate the cumulative effects over time. 
Some general impacts predicted by the integrated modelling system (Teck 
Resources, 2011) include: 


i. Drawdown propagation due to basal water 
depressurization will largely be constrained to the mining 
project area. Groundwater levels will begin to recover 
following shutdown of depressurization pumping, 
and the groundwater flow model predicts that far 
future groundwater levels would be similar to those of 
predevelopment. 


ii. Activities such as muskeg drainage and overburden 
dewatering during mine construction and operation 
will result in increased flows to receiving watercourses. 
Reductions in drainage area because of closed-circuit 
operation and the creation of pit lakes at closure will 
reduce flood flows to receiving waters. 


iii. Oil sands developments were predicted to have negligible 
effects on acute and chronic toxicity and tainting 
potential in all receiving waters in the local and regional 
study areas. The proposed mitigation measures will 
ensure that acute and chronic toxicity and tainting 
potential will be at levels appreciably lower than the 
corresponding guideline or threshold values, and that 
additional adaptive management options exist in the 
event that they are required. The concentrations of several 
substances are predicted to increase above base case 
snapshots but remain below guidelines or chronic effects 
benchmarks (CEBs). 
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iv. Fish habitat in the local study area (LSA) is primarily 
of low value and composed of forage fish typical of the 
Athabasca Oil Sands Region. Construction of oil sands 
projects will result in the alteration or destruction of 
fish habitat and an associated loss of fish relative to 
abundance, but it will have no effect on fish or fish habitat 
diversity. The loss of fish habitat will be compensated 
for by the construction of the lake, which will offset the 
potential loss of fish habitat. 


The ARM model, originally developed by Golder Associates and applied 
to EIAs of oil sands projects, is a two-dimensional vertically-averaged 
model based on an analytical solution to river dispersion equations under 
steady-state conditions, and implemented using VBA (Visual Basic for 
Applications) and the Microsoft Excel application. It has been updated by 
Four Elements Consulting Ltd. to include a new functionality for opti- 
mal regional substance load allocation in the Lower Athabasca River. 
Dynamic-link library (DLL) techniques were adopted to speed up com- 
puting time (Four Elements Consulting Ltd., 2014a, 2014b). A range of 
water quality parameters (including 11 general indicators such as chloride, 
TN and TP, 28 metals, 19 PAHs, total phenolics, and toxicity-chronic) can 
be simulated by the model. However, flow in the river was calculated using 
Leopold-Maddock equations rather than a hydrodynamic model. Similar 
to the DOSTOC model for the Lower Athabasca River, dynamic ice pro- 
cesses and river bed sediment transport processes cannot be simulated 
with the ARM configuration. 

More recently, an integrated hydrodynamic and water quality 
modelling framework for the Lower Athabasca River was proposed by 
Environment Canada and Climate Change (ECCC), consisting of MIKE 
11 for long-term one-dimensional simulations and EFDC for short-term 
detailed two-dimensional simulations, each externally coupled with a 
one-dimensional MIKE-ICE or CRISSP-1D model (Dibike et al., 2018; 
Kashyap et al., 2017; Shakibaeinia, Dibike et al., 2017; Shakibaeinia, 
Kashyap et al., 2016). A demonstrative rather than a full range of water 
quality parameters required for assessing impacts of oil sands development 
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are configured in the models, including TSS, BOD, DO, phosphorus, and 
nitrogen components, three PAHs (pyrene, phenanthrene, and Cl-benz[a] 
anthracenes/chrysenes) and three metals (lead, arsenic, and vanadium). 
The MIKE 11 model was developed for a 10-year historical period (2001- 
2010) and the EFDC model was calibrated over short periods under steady 
flow conditions. Both models were applied to hypothetical future scen- 
arios. Non-point source loadings were estimated from limited measure- 
ments, rather than from watershed modelling of land use changes and de- 
velopment scenarios. The modelling framework is under further develop- 
ment to enhance the configuration and functionalities by restructuring 
the model grid for EFDC and cross-section profiles for MIKE 11, based 
on derived high-resolution DEM (Chowdhury, 2017). A two-dimensional 
hydrodynamic, sediment transport and water quality model for the Lower 
Athabasca River has been developed and calibrated using EFDC+ with a 
number of enhancements, including: extended model domain to cover the 
main channel and 10-year floodplain from the Athabasca River upstream 
of Fort McMurray to the Athabasca River near Old Fort; optimized mod- 
el grid to improve computational burden to allow longer period (2000- 
2016) simulations in reasonable run-time; improved hydrodynamic model 
simulating dynamic ice cover formation and melting processes; cohesive 
and noncohesive sediment transport with the best available TSS and river- 
bed sediment data; and improved water quality (eutrophication and DO) 
and toxic (three representative toxics) modules by including all major 
tributaries, withdrawals and returns on the Lower Athabasca River, and 
point source effluents in the oil sands region (DSI, 2019). The river models 
are expected to be coupled with a distributed watershed model to achieve 
integrated watershed and river water quality modelling, as part of a multi- 
year mission to develop a comprehensive and integrated environmental 
modelling system for the cumulative effects assessment in oil sands region. 


5.0 Modelling in the Athabasca River Watershed - Case Study 129 


5.2 Atmospheric deposition and acidification 
modelling in the Athabasca Region 


In 2003, RWDI West Inc. (2003) conducted a CALMET/CALPUFF mod- 
elling study for CEMA. The model was run for a one-year period to pre- 
dict ambient concentrations and annual deposition rates for 39 priority 
substances in the Athabasca Oil Sands Region, that can then be used to 
screen potential human health risks using multimedia risk assessment 
techniques. The 39 priority substances, including SO,, NO,, and VOCs, 
were identified in the emission inventory of the oil sands region. Based 
on comparisons with available monitoring data, the highest level-of-con- 
fidence in model predictions is associated with the priority substances the 
emissions of which are well-defined, whereas a lower level-of-confidence 
is associated with emissions of priority substances from fugitive sources. 
For some contaminant-receptor combinations, the predicted ambient lev- 
els are so sufficiently low that a background term is required to allow for a 
representative comparison. 

CEMA considered applying one or both of two widely used air quality 
modelling systems, CALPUFF and/or CMAQ, for sulphur and nitrogen 
deposition modelling to assess historical, current, and future environ- 
mental exposures due to emissions from the oil sands industry and 
other sources in and around the Regional Municipality of Wood Buffalo 
(RMWB). ENVIRON International Corporation and Stantec Consulting 
Ltd. (2012) conducted a detailed evaluation and comparison of the per- 
formances of the CALPUFF and CMAQ modelling systems to under- 
stand the relative strengths and weaknesses of the two models for imple- 
mentation of the Acid Deposition Management Framework (ADMEF). It 
was determined that the differences between the two model outputs are 
strongly related to differences in the model inputs in spite of sharing the 
same data sources. There is no basis for choosing one model (CALPUFF 
vs. CMAQ) over the other. Both may be used, each with different advan- 
tages. The CALPUFF model has the advantage of consistency with previ- 
ous ADMF studies and other EIA studies, predicts the peak one-hour SO, 
and NO, concentrations better than CMAQ, and requires a lower level of 
effort to apply. The CMAQ model has the advantage of improved chem- 
istry algorithms and appears to perform better in simulating lower SO, 
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and NO, concentrations and the wet deposition of sulphur and nitrogen 
compounds. 

CALPUFF modelling (Exponent Inc., 2014) was also conducted to es- 
timate acid deposition to provide input to the Model of Acidification of 
Groundwater in Catchments (MAGIC) to determine whether acid depos- 
ition-related changes to lakes and soils will or will not exceed defined 
thresholds. The CALPUFF model was calibrated with 2010 meteorology 
and observation data and then applied with 1980 representative year me- 
teorology data, as well as with historical, base, and two future emission 
scenarios. The model was run for a one-year period to predict the annual 
dry, wet, and total depositions of acidifying sulphur and nitrogen com- 
pounds. The deposition measurement data was considered insufficient for 
a reliable comparison between measurements and predictions. In general, 
the predicted total Potential Acid Input (PAI) in 1980 is slightly lower than 
that in 2010 for the existing emission inventory. The area of peak total 
PAI is located in the Fort McKay area, which is the area with the largest 
emissions. The second highest peak area is located in the Bridge View area 
along the southern boundary of the domain. 

A CMAQ modelling study was completed with the same source of 
data and procedures as Exponent’s CALPUFF modelling (S. Cho et al., 
2017). Similar findings were obtained that include (i) a predicted area of 
high sulphur and nitrogen deposition near the largest oil sands operations 
in Alberta’s oil sands region, and predicted higher dry deposition than 
wet deposition in the study area; (ii) the predicted gross PAI (wet & dry) 
deposition increases from the historical to an existing case with further 
increases for the two future scenarios; and (iii) the nitrogen deposition 
predicted by the model, which comprises, on average, approximately 60% 
of the total PAI acidic deposition in the region. 

ECCC calculated estimates of potential effects on ecosystems in 
the Canadian provinces of Alberta and Saskatchewan due to acidify- 
ing deposition (Makar et al., 2018). Based on a one-year simulation of a 
high-resolution implementation of the Global Environmental Multiscale- 
Modelling Air-quality and Chemistry (GEM-MACH) model, the critical 
loads of sulphur and nitrogen deposition (dry, wet, and total) for aquatic 
and terrestrial ecosystems were derived. The spatial extent of the regions 
exceeding critical loads varied between 1E+4 and 3.3E+5 km’, for the 
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more conservative observation-corrected estimates of deposition, with 
the variation dependent on the ecosystem and the critical load calcula- 
tion methodology. Other findings of the study were outlined and include 
(i) the evaluation of the model simulation against two different sources 
of deposition data - total deposition in precipitation and total deposition 
to snowpack in the vicinity of the Athabasca oil sands, (ii) the variabil- 
ity of observed ions in wet deposition in precipitation (observed versus 
model sulphur, nitrogen, and base cation R? values of 0.90, 0.76 and 0.72, 
respectively), while being biased high for sulphur deposition, and low for 
nitrogen and base cations (slopes 2.2, 0.89 and 0.40, respectively), and (iii) 
the predicted potential ecosystem effects within each of the regions, repre- 
sented by the ecosystem critical load datasets using a combination of 2011 
and 2013 emissions inventories. 

The Model of Acidification of Groundwater in Catchments (MAGIC) 
has been applied to soils and lake catchments in the oil sands region to 
determine sensitivity to acid deposition under two deposition scenarios 
(base case and double acid) (Whitfield et al. 2009, 2010, 2011; Whitfield 
and Watmough, 2010). The model simulated average monthly or annual 
soil solution and surface water concentrations for sulfate (SO,”), calcium 
(Ca**), magnesium (Mg”*), sodium (Na*), potassium (K*) and pH, as well 
as exchangeable soil fractions of Ca**, Mg”*, Nat, and K*. Lumped indi- 
cators include base cation (BC) and acid neutralizing capacity (ANC) for 
lakes, base saturation (%) and critical threshold, and molar base cation 
to aluminum ratio (BC:Al) for soil physicochemical characterization. 
The research demonstrated that lakes in the region would not be at risk 
of acidification under either deposition scenario. In contrast, forest soil 
weathering rates range from very low to moderate, and soil chemistry was 
predicted to change under both deposition scenarios. 

With the annual atmospheric deposition rates generated by CALPUFF 
modelling, a modelling approach was implemented within Golder’s com- 
prehensive oil sands environmental modelling system for predicting the 
contributions and potential effects of aerially deposited PAHs and metals, 
as well as the impacts of snowpack and snowmelt on water quality from 
proposed oil sands developments (Dayyani, Daly, & Vandenberg, 2016). 
The contribution (loading) of snowmelt to surface water concentrations 
was estimated using two methods: a conservative mass balance approach 
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adopted for metals, and an unsteady-state, mass balance, and multi-com- 
partment fate model, or Coastal Zone Model for Persistent Organic 
Pollutants (CoZMo-POP). Simulations were run with a timestep of one 
hour. The model predicted that: i) under base case conditions, cadmium 
and chromium concentrations were both predicted to exceed guideline 
values in the Muskeg River; ii) All PAHs were predicted to remain below 
guidelines in surface water for both projects and under all assessment 
cases. 

The model was thought to over-predict the concentration of metals in 
receiving water pathways because it did not account for retention of met- 
als during snowmelt in the soil matrix. Refinements to the mass balance 
approach for modelling metals should focus on retention of metals on the 
landscape during the melt period. Refinements to CoZMo-POP might in- 
clude application of the water-sediment partitioning module to predict 
sediment PAH concentrations. 


5.3 Watershed modelling in the Athabasca 
Region 


A variety of watershed models have been used in the region such as HSPF, 
SWAT, VIC, WATFOOD, and MISBA. Pietroniro et al. (2006) investigated 
the potential effects of climate change on the hydrological regimes of three 
large lakes and two inflow sources (from the Peace and Athabasca riv- 
ers) on the Peace-Athabasca Delta. They indicated that changes in climate 
can result in an earlier melt season, higher winter flows, an increase of 
flows (more at headwater), and a reduction of peak flow. Kerkhoven and 
Gan (2006, 2011) assessed the impact of climate change using the MISBA 
model. They reported a large decline in spring snowpack, annual runoff 
(-21%), mean maximum annual flow (-4.4%), and mean minimum flow 
(-41%) by the end of century for ARB. The decline in the average annual 
flows has been also reported by Golder Associates (2009) when they evalu- 
ated the climate impacts using the HSPF model. An increase in winter 
flows and a decrease in summer flows have been found in Eum etal. (2014) 
and in Leong and Donner (2015), through assessing the impacts of cli- 
mate change on the hydrology of the region using VIC and Integrated 
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Biosphere Simulator - Terrestrial Hydrology Model with Biogeochemistry 
(IBIS-THMB), respectively. 


5.4 Groundwater modelling in the Athabasca 
Region 


CEA for groundwater requires a groundwater model with a mesh/grid 
that can conform to complex geology and hydrogeology and with the cap- 
ability of simulating density-dependent flow and transport, as well as sur- 
face water-groundwater interactions. A group of 10 in situ oil sand oper- 
ators (the SAOS Group) initiated a process in 2007 to develop a regional 
surface water-groundwater model for the oil sands in situ area south of 
Fort McMurray. They conducted a review of surface water-groundwater 
numerical modelling to select the most suitable model. MODFLOW and 
FEFLOW were selected with the capability of being coupled with sur- 
face water models — as MODHMS and FEFLOW/MIKE 11, respectively 
(WorleyParsons, 2010). Major groundwater modelling studies have been 
undertaken in the Lower Athabasca River watershed. For instance, in 
2009, WorleyParsons Canada Services Ltd. developed a 3D MODFLOW 
groundwater model (2 km grid resolution) for the entire Lower Athabasca 
Regional Plan (LARP) region to assess the potential impacts from oil 
sands development. They indicated that the current drawdown for the 
various major aquifers was substantial and that the trend of drawdown 
was anticipated to decrease. Subsequently, modelling studies have been 
more regionally focused and based on operational purposes whereas the 
Athabasca Oil Sands (AOS) area has been split into the following three 
regions: 


1. Northern Athabasca Oil Sands (NAOS) region, which 
mainly contains the surface mineable deposits. The 
boundaries of the NAOS study area are 


e North: starting in the northwest, the boundary follows 
the Sand River, Athabasca River, Firebag River, and 
Marguerite River sub-basin; 

e East: Alberta/Saskatchewan border; 
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e South: Athabasca and Clearwater Rivers; and 

e West: northwest boundary follows the western extents 
of the Gardiner Lake and Snipe Creek sub-basins, 
then follows the Dunkirk River south to the western 
boundary of the MacKay River sub-basin, south to the 
Athabasca River. 


2. Southern (SAOS) area where in situ extraction occurs in 
the region. The boundaries of the SAOS region are: 


e North: Athabasca and Clearwater Rivers; 

e East: Alberta/Saskatchewan border and Christina River 
Sub-basin; 

e South: Centre of Township (T) 69 from Range (R) 1 to 9 
West of the Fourth Meridian (W4M) continuing along 
the Beaver River Basin; and 

e West: Southwest boundary follows the La Biche sub- 
basin to the confluence of the Athabasca and La 
Biche River. The boundary continues north along the 
Athabasca River. 


3. Cold Lake region which has been given the name Cold 
Lake-Beaver River (CLBR) based on its location within 
the Beaver River Basin. 


For the NAOS region, WorleyParsons Canada Services Ltd. developed a 
3D FEFLOW model in 2012, which allows for a flexible mesh refinement, 
simulates density-dependent and fractured flow, and includes river flow 
interactions by linking the model to MIKE 11. The first phase of the model 
development was mainly focused on the configuration of the model and 
understanding the hydrology/hydrogeology of the region. They reported 
that the recharge rate ranges from 247 to 1,150 million m*/year for the 
NAOS study area. The total annual discharge from groundwater to the 
Athabasca River (including its tributaries) ranges from 236 million m°/ 
year to 590 million m*/year. In addition, as of October 2011, it was esti- 
mated that the rate of groundwater withdrawal was roughly 41.5 million 
m’/year within the NAOS region. A FEFLOW model has been also used 
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by WorleyParsons to model groundwater of the SAOS region in 2010. They 
reported that the recharge in the SAOS ranges from 290 million m*/year 
to 1,500 million m?/year. The total discharge from groundwater to the 
river ranges from 540 million to 1,350 million m*/year. As of February 
2009, it was estimated that the total annual non-saline allocation volume 
is more than 15 million m° within the SAOS region. In 2016, the FEFLOW 
model developed by WorleyParsons for the SAOS was further modified by 
MATRIX Solutions Inc., to improve the calibration of the model. 

For the CLBR region, a MODFLOW groundwater model was de- 
veloped by the Alberta Geological Survey in 2005 to better understand the 
regional water balances and groundwater flow regimes. It was reported 
that recharge rates estimated by the model were highest in the northeast 
and southeast of the domain area (13.2 mm/year and 7.6 mm/year, re- 
spectively) and were lower than 5 mm/year in the western portion. 


5.5 Surface water and groundwater interactions 
in the Athabasca Region 


Few attempts have been made to model surface water-groundwater (SW- 
GW) interactions. In a study conducted by WorleyParsons (Integrated 
Sustainability Consultants Ltd., 2013), the MODFLOW model was used 
for the lower Athabasca regional planning region and they found 50% 
exceedance of available drawdown in the basal McMurray Formation, 
especially in the mineable area. They also indicated that a more compre- 
hensive modelling tool is required for CEA under multiple-stressors. A 
review of modelling studies for assessing the potential cumulative impacts 
to groundwater and surface water in the MacKay River watershed was 
conducted by the Cumulative Environmental Management Association 
(CEMA, 2014). It was predicted that groundwater discharge can reduce 
to -0.001 to 0.043 m*/s when various oil sands projects are developed 
compared to the current condition 0.01-0.055 m?/s. They indicated that 
many project EIAs used different models with different data sets and as- 
sumptions, which makes comparison of projects and impacts difficult. 
Kassenaar (2016) assessed the potential cumulative impacts to SW-GW 
from in-situ oil sands operations using GSFLOW in the MacKay River 
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watershed. They indicated that drawdowns do not (on a watershed-scale) 
appear to grow over time. However, cumulative groundwater diversions 
appeared to create unsustainable local impacts under extreme and defined 
scenarios. Furthermore, the simulations showed that that groundwater di- 
versions may significantly affect small to intermediate sized tributaries. 


5.6 Land use/land cover modelling in the 
Athabasca Region 


A State-and-Transition Model (STM) model was developed for CEMA 
and validated by Apex Resource Management Solutions Ltd. to simulate 
reclamation dynamics for the mineable oil sands region of Alberta, such 
that it supports the development of a Reclamation Classification System 
(RCS) (Daniel, 2011; Frid & Daniel, 2012). The model utilized landscape 
vegetation state and the probability of transition between states to find 
potential probabilities of changes. STMs can help managers establish a 
land classification system by describing the land units and phases in the 
system, including the relationships between the various units and phas- 
es. STMs explicitly recognize the relationship between management al- 
ternatives and land classification, ensuring that the classification system is 
management-oriented. A quantitative STM was developed which can be 
used to both conceptually and quantitatively model landscape level chan- 
ges over time because of alternative reclamation scenarios. The model was 
parameterized with existing data from the Long-Term Plot Network and 
other CEMA projects. 

The Alberta Biodiversity Monitoring Institute (ABMI) is one of the 
main sources of land use and land cover data which also provides new 
analytical methods and visualization approaches to deliver geospatial 
products, such as the ABMI province-wide wall-to-wall Human Footprint 
Inventory (HFI). The ABMI has also developed a predictive product, 
called Predictive Landcover (PLC), which classifies land cover into seven 
classes (open water, bog, marsh, swamp, fen, upland, and wetland general). 

Another important product to support land use and land cover deci- 
sion-making in Athabasca is ALCES, a landscape and mapping software. 
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ALCES has been applied to inform planners and stakeholders about pos- 
sible future outcomes associated with land use and development. 


5.7 Climate change in the Athabasca Region 


Climate change can be considered as a natural stressor which needs to be 
taken into account along with other stressors for CEA. Climate change 
studies in the Athabasca watershed can be divided into two main categor- 
ies: the assessment of climate change and variability, and the impact of 
climate change on environmental processes. Regarding climate change 
and variability in the Athabasca watershed, Bonsal and Cuell (2017) indi- 
cated that periodic extreme droughts and excessive moisture conditions 
are expected mainly due to persistent and mid-tropospheric circulation 
patterns that disrupt expected temperature and precipitation in the region 
(Bonsal & Cuell, 2017). Furthermore, substantial inter-annual variability 
with more drought-like summer and slightly wetter annual conditions, as 
well as decadal-scale variability are predicted over the entire region. On 
the other hand, according to the climate change impact studies, a com- 
parison between current and future climate conditions in the watershed 
has indicated a significant shift towards an earlier melt season and also an 
increase in winter flows (Toth et al., 2006). The shortened snowfall sea- 
son along with an increase in sublimation can result in a decline in the 
spring snowpack as well as an expected decline in average annual flows 
(Kerkhoven & Gan, 2011). 


5.8 Limitation of modelling for cumulative effects 
assessment in the Athabasca Region 


In order to build an understanding of changes in the environment and 
to ensure that the cumulative effects of the stressor/s are captured, CEA 
often requires to take into account: (i) the related theme and sub-theme 
areas and their environmental processes, (ii) interactions and feedback 
mechanisms between environmental processes (in different related theme 
and sub-theme areas), and (iii) both local and regional response scales. 
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However, the majority of modelling studies in Athabasca have the follow- 
ing limitations: 


1. They often take into account either limited theme areas or 
environmental processes. 


2. ‘They often lack the capacity to capture the interactions 
between major system components. 


3. They only consider either regional-scale response (which 
lacks the necessary detail on local-scale pathways) or 
local response (which might not be suitable for guiding a 
regional-scale decision-making process). 
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6.0 INTEGRATED MODELLING 
FRAMEWORK FOR CEA 


This book has presented a comprehensive review of literature on environ- 
mental components, how they are connected and interact with each other, 
and the application of numerous modelling tools to understand the com- 
plex behaviors of environmental factors. Where possible, shortcomings of 
model applications have been highlighted to better understand the limita- 
tions and critical gaps in modelling the environment in a holistic way. The 
integrated environmental modelling framework proposed here addresses 
the gaps and limitations identified through this comprehensive literature 
review and could be applied at various temporal and spatial scales. The 
developed framework consists of a core and three supported layers, de- 
scribed as follows (Figure 9). 


Core: The core consists of a comprehensive integrated watershed model- 
ling system with the following characteristics: 


1. a physically based model which provides a detailed 
description of the processes that occur in the watershed; 


2. a fully distributed model which considers the watershed 
as finite geo-referenced computational units with 
different responses to forced inputs; 


3. adynamic model that can employ both a short- and long- 
timestep along with either a detailed or coarse drainage 
network schematization of the watershed; 
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4. adynamic model that can incorporate land surface 
changes at different time intervals during the simulation 
of environmental processes; 


5. amodel that can simulate environmental processes for 
both local and regional scales; 


6. a model which either includes all the required processes 
to assess water (surface and subsurface) quantity, quality, 
land, climate, ecology, and air deposition components, 
or can be linked to other models to further evaluate the 
required mentioned processes; and 


7. amodel that can be linked to a socio-hydrology, -ecology, 
-economic model, as well as a decision support system. 


Layer 1: Although, the core modelling system is capable of simulating sur- 
face water-groundwater at local to regional scale, this layer is considered 
when further assessment - with a higher resolution - is required for either 
surface water-groundwater or ecological components. This can support 
the modelling system for specific applications such as when strong in- 
homogeneities are present in the modelling domain or when a detailed 
cumulative effects assessment is required which is also linked to a regional 
CEA assessment. 


Layer 2: this layer of the proposed integrated environmental modelling 
framework supports predictive modelling through scenario building for 
various system drivers, such as climate change, land use, atmospheric 
depositions, management practices, or policy changes. Scenarios de- 
veloped under this layer are simulated using the modelling core and/or 
layer 1 to develop future or what-if projections. 
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Figure 9. Integrated Mechanistic Based Modelling Framework for CEA 


Layer 3: this layer incorporates the interrelations between the social and 
economic factors with the environment to fill the gaps between policy 
development (public policy, economic policy, environmental policy, etc), 
decision-making, and science. This layer supports the modelling system 
to capture the interactions between three sustainability pillars (i.e. eco- 
nomic, environmental, and social) essential for existence of humankind 
on planet earth. The interactions could ideally be through a web-based 
interface which offers the flexibility of investigating a wide range of scen- 
arios relating people and economy to environment. 
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6.1 Coupling strategy 


A fully dynamic coupling is desired for the core of the modelling system 
which would dynamically capture the interactions and exchanges occur- 
ring with sub-systems of environment (i.e., land, surface water, ground 
water, air). Models on different layers (1, 2, and 3) could be linked using 
appropriate coupling strategy e.g., dynamic coupling (e.g., receiving water- 
body hydrodynamic model and water quality model), iterative coupling 
(e.g., optimal water and load allocation model and receiving water quality 
model), or interactive coupling (e.g., land use model for changing time 
periods and watershed model). The linkages between core and other layers 
of the proposed framework could be achieved through hybrid coupling 
strategy where most models in different domains or disciplines are inte- 
grated through a sequential coupling. 


6.2 Selection of models 


The models should be selected based on their ability to capture the cumu- 
lative effects of identified stresses on defined environmental processes. 
Their strengths in respective disciplines should be considered as well as 
their acceptability as the state-of-the-art model with demonstrated ap- 
plications worldwide. Furthermore, factors other than model capabilities 
should also be considered, such as the desired degree of complexity for 
cumulative effects assessment, watershed characteristics, and temporal 
and spatial scale. 

The capability of the integrated modelling system is highly reliant 
on the selected core system and its characteristics (see section 6.1). Based 
on our literature review, the five most popular models, which might be 
able to meet the core requirements, are: GSFLOW, MIKE SHE/MIKE 
11, HydroGeoSphere, MODHMS, and ParFlow. An appropriate model or 
models should be selected based on their capabilities and the objectives of 
CEA projects. 
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6.3 Novelty of the proposed integrated 
environmental modelling framework 


Despite the conventional cumulative effects modelling approach which 
simplifies representation of environmental processes within a single en- 
vironmental media, the proposed integrated environmental modelling 
framework brings together a set of interdependent components by char- 
acterizing the stress—response relationships based on modelling of inter- 
actions of a variety of components and cross-pathways. In other words, 
the proposed framework introduces a holistic systems-based approach 
that integrates multidisciplinary environmental components that can fa- 
cilitate cumulative effects management strategies. The framework has the 
following main charactristics: 


e It can explore dynamic, nonlinear, and complex 
interactions and feedbacks among environmental 
processes. 


e It can tease out impacts of multiple stressors on 
environmental processes. 


e It is capable of simulating environmental processes at 
different spatial and temporal scales. 


e It serves the cumulative effects management needs 
to understand the complexity of the system by 
involving stakeholders and scenarios development, and 
consequently analyzing trade-offs among alternatives. 


6.4 The challenges 


CEA integrates various meteorological, hydrological, hydro-geological, 
biological, and chemical processes that occur at a broad range of spatial 
and temporal scales that characterize environmental processes. For in- 
stance, infiltration and soil hydraulic conductivity are pore scale and field 
scale processes, respectively, while transpiration occurs at a leaf area scale. 
On the other hand, overland flow is a watershed (local/regional) scale pro- 
cess, while climate change processes occur at the global scale. In terms 
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of temporal scale, turbulent fluxes (e.g., moisture, sensible heat) and the 
land-atmosphere exchanges occur at the scale of seconds to minutes, while 
changes in land surfaces occur on the order of years to decades. Although 
interactions between these processes at different scales are well-estab- 
lished observationally and theoretically in integrated physical modelling 
systems, computational limitations have restricted the use of integrated 
cumulative effects predictive models for local and regional studies. For 
example, there are still many challenges, such as adaptation of water re- 
sources to climate and anthropogenic stressors that need to be addressed 
across large domain scales with a fine resolution, but integrated cumula- 
tive effects predictive modeling for a large domain scale is an intractable 
computational problem. 
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Investigating the complex nature of environmental problems requires the 
integration of different environmental processes across major components 
of the environment. Cumulative effects assessment (CEA) not only includes 
analyzing and modelling environmental changes, but also supports planning 
alternatives that promote environmental monitoring and management. 


The adoption of integrated modelling approaches requires the development 
of frameworks which may be used to investigate individual environmental 
processes and their interactions with each other. Integrated modelling 
frameworks are often the only way to examine important environmental 
processes and interactions, relevant spatial and temporal scales, and feedback 
mechanisms of complex systems for CEA. 


This book examines the ways in which interactions and relationships between 
environmental components are understood, paying special attention to climate, 
land, water quantity and quality, and both anthropogenic and natural stressors. 
It reviews modelling approaches for each component and existing integrated 
modelling systems for CEA. Finally, it proposes an integrated modelling 
framework and provides perspectives on future research avenues for 
cumulative effects assessment. 
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